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SECTION  1 

NUMERICAL  SIMULATION  OF  2-D  FLOWS  USING  A  2-D 
ADAPTIVE  FINITE  ELEMENT  SCHEME 


A  transient,  two-dimensional,  finite-element  shock-capturing  scheme  on  unstructured 
grids  was  applied  to  the  study  of  a  shock  dif&acting  siround  a  2-D  structure,  representing 
a  2-D  cut  of  a  3-D  train,  suspended  above  the  rails  and  a  rigid  elevated  surface.  The  area 
between  the  trmn  and  the  surface  was  partially  blocked  by  the  rails,  resulting  in  com¬ 
plex  shock  diffraction  processes.  The  resiilts  demonstrate  the  capability  of  the  developed 
adaptive  refinement /coarsening  algorithm  to  properly  adapt  to  weak  shocks,  expansions 
and  contact  discontinuities,  and  highlight  the  resulting  excellent  resolution  of  the  captured 
flow  features.  In  addition  to  interesting  shock  diffraction  and  propagation  phenomena,  the 
results  demonstrate  the  capability  of  the  new  code  to  capture,  and  define  in  great  detail, 
vortex  sheets  shed  from  sharp  comers.  We  show  that  the  baroclinic  effect,  an  inviscid  pro¬ 
cess,  controls  the  shedding  phenomenon  during  the  diffraction  phase.  Hence,  the  Eulerian 
model  is  able  to  correctly  predict  this  process. 

1,1  INTRODUCTION. 

During  the  past  ten  years  the  CFD  community  has  experienced  a  proliferation  of 
shock  capturing  schemes  whose  ultimate  objective  is  the  sharp,  nonoscillatory  capturing 
of  transient  shocks,  even  those  that  have  propagated  great  distances  [1*3].  The  practical 
objective  of  these  methodologies  is  to  predict  the  loads  exerted  by  shocks,  initiated  by  nu¬ 
clear  or  High- Explosive  binsts,  on  stationary  or  moving  structures  located  a  long  distance 
from  biust  point.  The  results  of  such  simulations  should  help  improve  targeting  specifica¬ 
tion  requirements  (for  offensive  purposes)  or  hardening  criteria  (for  defensive  objectives). 
FVom  the  numerical  algorithm  development  point  of  view,  the  long  propagation  times  and 
distances  that  the  shocks  are  likely  to  transverse  from  burst  to  impact  p>oint,  j>ose  signif¬ 
icantly  greater  demands  than  required  from  traditional  shock  capturing  schemes,  which 
are  only  intended  to  produce  steady,  sharp,  nonoscillatory  shocks  on  airfoils,  etc.  [4-6] 
to  }deld  the  steady  lift  and  drag  forces.  Past  simulations  of  shock  wave  propagation  over 
long  distances  relied  on  either  fixed-mesh  structured  grid  approach  [2]  or  sliding  refined 
zones,  intended  to  continously  surround  the  shocks  with  a  finer  mesh  than  elsewhere  in 
the  domain  [7].  The  first  approach  yielded  coarser,  diffused  shocks  after  some  propaga¬ 
tion  distance,  since  even  for  two-dimensional  calculations  it  is  not  economically  feasible, 
even  with  present  day  class  VI  computers,  to  grid  a  domain  of  several  himdreds  or  thou¬ 
sands  of  meters  with  1  cm  size  cells.  The  second  approach  was  much  more  economical, 
as  only  part  of  the  domain  was  finely  refined.  Nevertheless,  the  computational  resources 
required  for  two-dimensional  computations  (both  memory  and  computational  time)  were 


1 


very  high,  while  three-dimensional  computations  would  have  been  prohibitively  expensive 
for  production- type  nms.  Hence,  the  objective  of  the  present  effort  was  a  natural  evolution 
of  past  deficiencies,  i.e.,  develop  a  numerical  scheme  capable  of  dynamically  adapting  to 
traveling  shocks  and  other  flow  discontinuities.  This  approach  allows  us  to  use  coarse  grid 
resolution  everywhere  in  the  computational  domain  where  flow  gradients  are  low,  while  ap¬ 
plying  very  fine  grid  resolution  wherever  better  flow  resolution  is  required.  Hence,  shocks 
should  always  advance  into  finely  refined  zones,  while  grid  coarsening  should  be  obtsiined 
in  areas  already  traversed  by  the  shocks. 

The  section  will  first  describe  the  CFD  Methodology  employed,  starting  with  the  de¬ 
sign  criteria  used  for  the  flow-solver.  Thereafter,  the  adaptive  refinement  scheme  employed 
is  presented,  followed  by  an  extensive  description  of  the  results  obtained  for  the  shock-train 
case  studied.  The  accent  of  the  section  is  on  the  results  obtained  imder  the  project  rather 
than  on  the  algorithms  employed  (which  were  not  developed  under  this  effort).  Therefore, 
the  expK>sition  of  the  algorithms  is  kept  to  a  sufi^cient,  but  not  exhaustive,  depth. 

1.2  CFD  METHODOLOGY. 

The  objective  of  the  present  effort  was  to  simulate  many  2-D  and  3-D  shock  diffraction 
processes  aroimd  complex-geometry  structures.  Thus,  the  flow  solver  employed  must: 

1)  Be  based  on  grid  systems  that  can  discretize  domains  of  arbitrary  complexity; 

2)  Be  able  to  simulate  moving  or  stationary  shocks  without  spurious  overshoots; 

3)  Be  able  to  allocate  in  the  most  efficient  maimer  the  degrees  of  freedom  or  equivalent 
computational  resources  during  the  course  of  a  simulation.  In  the  present  case,  these  three 
design  criteria  are  met  as  follows: 

1)  In  order  to  discretize  domains  of  arbitrary  complexity,  unstructured  grids  in  conjrmction 
with  Finite-Element  Methods  are  used. 

2)  In  order  to  simulate  accurately  the  strong  shocks  present  in  the  flowfields,  a  high- 
resolution  monotonicity-preserving  algorithm  for  unstructured  grids  is  used.  The  method, 
called  FEM-FCT,  is  based  on  Zalesak’s  [8]  generalization  of  the  Flux- Corrected  Transport 
(FCT)  algorithms  of  Boris  and  Book  [9,10]  to  multidimensional  problems. 

3)  In  order  to  allocate  in  the  most  efficient  manner  the  degrees  of  freedom  of  the  mesh,  an 
efficient  adaptive  refinement  technique  for  transient  problems  is  employed. 

1.2.1  The  Flow  Solver;  FEM-FCT. 

As  stated  above,  high-resolution  monotonicity-preserving  schemes  must  be  developed 
to  simulate  the  strong  nonlinear  discontinuities  present  in  the  flows  under  consideration.  A 
number  of  these  schemes  have  been  develop>ed  over  the  last  years  [4-6,11,12].  The  scheme 
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used  here  is  based  on  Zalesak’s  generalization  (8]  of  the  1-D  FCT  schemes  of  Boris  and 
Book  [9,10].  Parrot  and  Christie  [13]  first  analyzed  FCT  schemes  in  the  context  of  Finite 
Element  Methods,  and  Lohner  et  al.  [14,15]  extended  these  ideas  further  to  include  the 
solution  of  systems  of  equations  and  the  consistent-mass  matrix  that  yields  high  temporal 
accuracy. 

Consider  a  set  of  conservation  laws  given  by  a  system  of  partial  differential  equations 
of  the  form 


dU  dFi  dFi 
dt  *  dx'  ~  dx‘  ’ 


(11) 


where  the  advective  fluxes  F«  =  Fa(U)  dominate  the  viscous  fluxes  F»  =  F^iU).  For  flows 
described  by  Eq.  1.1,  discontinuities  in  the  variables  may  arise  (e.g.,  shocks  or  contact 
discontinuities).  Any  numerical  scheme  of  order  higher  than  unity  will  produce  overshoots 
or  ripples  at  such  discontinuities  (the  so-called  Godunov  theorem).  Very  often,  particularly 
for  mildly  nonlinear  systems,  these  overshoots  can  be  tolerated.  However,  for  the  class  of 
problems  studied  here,  overshoots  will  eventually  lead  to  numerical  instability,  and  will 
therefore  have  to  be  suppressed. 

The  FCT  algorithm  combines  a  high-order  scheme  with  a  low-order  scheme  in  such 
a  way  that  in  regions  where  the  variables  under  consideration  vary  smoothly  (so  that  a 
Taylor  expansion  makes  sense)  the  high-order  scheme  is  used,  whereas  in  those  regions 
where  the  variables  vary  abruptly,  the  low-order  scheme  is  favored. 

The  temporal  discretization  of  E)q.  1.1  yields 


Ijn+I  (12) 

where  AU  is  the  increment  of  the  unknowns  obtained  for  a  given  scheme  at  time  t  =  t'*. 
Our  aim  is  to  obtain  a  AU  of  as  high  an  order  as  possible  without  introducing  overshoots. 
To  this  end,  we  rewrite  Eq.  1.2  as 

un+i  =  17"  +  Ai;'  -I-  {AU'*  -  AU’),  (1.3) 


or 


t7"+*  =  U’  +  {AU’'  -  AU’).  (1.4) 

Here  AU’'  and  AU’  denote  the  increments  obtained  by  some  high-  or  low-order  scheme, 
whereas  U’  is  the  (ripple-free)  solution  at  time  t  =  of  the  low-order  scheme.  The  idea 
behind  FCT  is  to  limit  the  second  term  on  the  right  side  of  Eq.  1.4; 
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C7"+*  =U*  +  Um{AU^  -  AU'), 


(1.5) 


in  such  a  way  that  no  new  overshoots  or  imdershoots  are  created.  In  addition,  strict 
conservation  on  the  discrete  level  must  be  maintained.  The  simplest  way  to  guarantee  strict 
conservation  for  node-centered  schemes  (only  these  are  considered  here)  is  by  constructing 
schemes  for  which  the  sum  of  the  contributions  of  each  individual  element  (cell)  to  its 
surrounding  nodes  vanishes  (“all  that  comes  in  goes  out”).  This  means  that  the  limiting 
process  (Elq.  1.5)  will  have  to  be  carried  out  in  the  elements  (cells).  FWther  details  may 
be  found  in  References  14  and  15. 

1.2.2  Transient  Adaptive  Refinement. 

A  veiy  attractive  featture  of  schemes  based  on  unstructured  grids  is  the  ease  with 
which  they  incorporate  adaptive  refinement.  The  addition  of  further  degrees  of  freedom 
does  not  destroy  any  previous  structure.  Thus,  the  flow  solver  requires  no  further  modifi¬ 
cations  when  operating  on  an  adapted  grid.  For  many  practical  problems,  the  regions  that 
need  to  be  refined  are  extremely  small  compared  with  the  overall  computational  domain. 
Therefore,  storage  and  CPU  requirements  are  typically  reduced  by  a  factor  of  10  to  100 
when  compared  to  an  overall  fine-resolution  fixed  mesh  [16-22].  It  has  been  our  experience 
that,  for  production  runs,  adaptive  refinement  has  been  the  crucial  element  necessary  to 
perform  simulations  at  an  acceptable  level  of  accuracy  in  a  reasonable  amount  of  time. 
The  simulation  presented  in  this  section  is  a  typical  example  of  such  runs. 

Several  authors  have  studied  adaptive  refinement  schemes  (7,16,23-25).  When  devel¬ 
oping  an  efficient  adaptation  methodology  for  transient  problems,  further  constraints  have 
to  be  considered:  a)  The  method  must  be  conservative,  i.e.,  a  mesh  change  should  not 
result  in  the  production  or  loss  of  mass,  momentum  or  energy,  as  this  will  produce  er¬ 
roneous  shock  propagation  speed  or  will  diffuse  shocks  or  contact-discontinuities;  b)  The 
method  should  not  produce  elements  that  are  too  small,  as  this  will  significantly  reduce 
the  allowable  time  step  of  the  explicit  flow  solver;  c)  The  method  must  be  fast,  as  the  grid 
is  modified  many  times  during  the  computation.  In  particular,  it  should  lend  itself  to  some 
degree  of  parallelism;  and  d)  The  method  must  not  be  memory  (storage)  intensive. 

These  constraints  are  met  by:  a)  Applying  strict  conservation  on  the  discrete  level  and 
conducting  the  limiting  process  of  FCT  on  the  elements;  b)  Using  classical  h-enrichment/ 
coarsening,  as  it  does  not  require  a  major  storage  overhead  and,  due  to  its  simplicity, 
lends  itself  easily  to  vectorization;  c)  Allowing  only  one  level  of  refinement/coar  ning  per 
mesh  change  in  order  to  minimize  the  logic  involved  and  thris  CPU  requirements,  and  d) 
Avoiding  successive  subdivision  of  a  triangle  into  two. 
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1.2.3  Grid  Logic. 


When  constructing  the  algorithm  to  refine  or  coarsen  the  grid,  one  faces  the  usual 
decision  of  speed  versus  storage.  The  more  information  from  the  previous  grid  is  stored, 
the  faster  the  new  grid  may  be  constructed.  As  storage  requirement  minimization  was 
one  of  the  goals  of  the  present  research,  we  tried  to  keep  only  the  essential  information 
needed  between  mesh  changes  without  sacrificing  an  excessive  amount  of  CPU  time.  In 
the  present  case,  six  integer  locations  per  element  are  kept  in  order  to  identify  the  “parent” 
and  “son”  elements  of  any  element. 

The  first  three  integers  store  the  new  three  ndghbor  elements  (“sons”)  of  an  element 
that  has  been  subdivided  into  four  (the  center  element  of  the  four  is  kept  as  “parent”).  If 
the  element  has  been  subdivided  into  two  elements,  the  neighbor  element  is  stored  in  the 
first  integer,  whereas  the  remaining  second  and  third  integer  locations  for  this  element  are 
set  to  zero. 

The  element  from  which  the  present  element  originated  (the  parent  element)  is  stored 
in  the  fourth  integer.  If  the  parent  element  has  been  subdivided  into  two,  the  negative  par¬ 
ent  element  number  is  stored,  allowing  the  distinction  between  the  1:2  and  1:4  refinement 
cases. 

The  fifth  integer  denotes  the  local  side  of  the  parent  element  from  which  this  element 
came. 

Finally,  in  the  sixth  integer  location  the  refinement  level  is  remembered.  These  six 
integer  locations  per  element  are  sufficient  to  construct  further  refinements  or  to  reconstruct 
the  original  grid. 

1.2.4  Error  Indicator. 

Many  possible  error  indicators  have  been  suggested  in  the  literature  [15-19,  23-25]. 
Numerical  experience  indicates  that  all  perform  similarly.  However,  the  following  require¬ 
ments  must  be  met  for  the  present  application:  a)  The  error  indicator  must  be  fast;  b) 
The  error  indicator  should  be  dimensionless,  so  that  more  than  one  'key-variable'  [19] 
can  be  monitored  simultaneously;  c)  To  be  applicable  to  a  large  class  of  problems,  the 
error  estimator  should  be  boimded  (independent  of  the  solution),  so  that  preset  refine¬ 
ment/coarsening  tolerances  can  be  employed;  d)  As  the  feature  may  move  only  very  slowly 
or  come  to  a  standstill  (e.g.,  a  shock  entering  a  very  dense  region),  the  error  indicator  must 
also  be  reliable  for  steady  state  applications;  and  e)  The  error  indicator  should  mark  for 
refinement  not  only  regions  with  strong  shocks,  but  also  regions  of  weak  shocks,  contact 
discontinuities  and  other  weak  features  in  the  flow. 

A  classic  interpolation  estimates  [19]  used  for  steady  state  computations  [20-22]  has 
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been  modified  to  meet  these  requirements.  These  estimators  make  use  of  an  appropriate 
seminoim  for  the  detection  of  those  regions  needing  further  refinement  or  coarsening,  e.g., 
the  H2*seminorm  [16,20-22,24,25] 


llu-u*llo<c-^^l«l2^  (1*6) 

where  u  denotes  the  exact  and  the  approximate  solution,  c  is  a  mesh-size  independent 
constant,  h  is  the  characteristic  mesh  size,  and 


Second  derivatives  are  justified  here  because  the  shape  functions  used  in  the  finite  element 
discretization  are  linear.  Numerically,  the  second  derivatives  at  the  nodes  are  evaluated  first 
via  a  variational  statement  the  integral  (7)  is  apprcodmated  conservatively  by  evaluating  for 
each  element  the  maximum  second  derivative  at  the  associated  nodes.  For  linear  elements 
of  constant  length  h  in  1-D,  one  obtains  for  the  first  step  at  the  nodes: 

ej  =  h-^-\Ui+i-2^Ui  +  Ui.i\.  (1.8) 

This  error  indicator  is  dimensional,  and  thus  is  not  bounded  a  priori.  A  major  defect 
noted  in  numerical  experiments  was  that  this  error  indicator  tends  to  favor  strong  shocks. 
This  results  from  the  fact  that  second  derivatives  of  a  key  variable  u  were  utilized.  Since 
the  objective  of  the  present  study  is  to  develop  a  methodology  that  accurately  simulates 
weak  shocks  and  contact  discontiniiities,  the  interpolation  theory  error  indicator  given  by 
Eq.  (1.8)  was  modified  as  follows: 


Ei  = 


\Uj^i-2-Ui  +  Ui.i\ 

\Ui+i-Ut\  +  \Ui-Uj-i\  +  0 


(1.9) 


where 


«  =  f[|P,+i|+2-|t7,|  +  |l^,-.|l. 

Dividing  the  second  derivatives  by  the  absolute  value  of  the  first  derivatives  yields 
an  error  indicator  that  is  boimded  (0  <  £/  <  1),  dimensionless,  and  prevents  the  eating 
up  effect  of  strong  shocks.  In  addition,  the  non-dimensional  error  indicator  permits  the 
use  of  several  critical  parameters  simultaneously.  Thus,  it  is  possible  to  use  both  density 
and  vorticity  for  problems  requiring  the  simultaneous  tracking  of  shocks  and  vortices.  The 
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terms  following  e  in  1.9  are  added  as  a  noise  filter  in  order  not  to  refine  wiggles  or 
ripples  that  appear  due  to  loss  of  monotonicity.  The  value  for  e  thus  depends  on  the 
algorithm  chosen  to  solve  the  PDEs  that  model  the  physical  process  investigated.  The 
multidimensional  form  of  this  error  indicator  is  given  by 


= 


i:i,AlaN!,Njdavjy 


1^.41  ‘  <^y‘ ' 


(1.10) 


where  denotes  the  shape-function  of  node  I. 

After  the  values  of  the  error  indicators  in  the  elements  are  determined,  all  elements 
lying  above  a  preset  threshold  value  are  refined,  while  all  elements  lying  below  a  preset 
threshold  value  are  coarsened. 


1.3  RESULTS. 

The  transient,  two-dimensional,  finite-element  shock-capturing  scheme  on  unstruc¬ 
tured  grids  described  above  was  applied  in  this  research  effort  to  the  study  of  a  shock 
interacting  with  a  train  suspended  above  a  rigid  elevated  surface.  The  area  between  the 
train  and  the  surface  was  partially  blocked  by  the  rails. 

Figures  1-1  through  1-11  show  combinations  of  computational  mesh,  grid  refinement 
levels,  pressure,  density,  vorticity  and  Mach  nrunber  contours  at  several  times.  The  grid 
was  adapted  every  seven  time  steps.  Density  was  chosen  as  the  key  variable  for  the  error 
indicator.  The  incident  weak  shock  travels  from  left  to  right,  with  a  shock  Mach  nmnber 
of  1.3.  The  color  contour  plots  presented  used  256  contoius,  with  blue  representing  the 
lowest  value  and  magenta  the  highest.  With  respect  to  the  mesh  refinement  procedure, 
one  level  of  mesh  refinement  is  defined  as  the  division  of  a  triangle  into  four,  performed  by 
connecting  the  mid-faces.  The  color  designations  for  the  mesh  refinement  levels  are  blue 
for  the  original  grid,  and  green,  yellow,  red,  magenta,  and  cyan  for  one,  two,  three,  foiu, 
and  five  refinement  levels,  respectively.  Two  parameters  control  mesh  refinement:  a)  the 
maximum  number  of  refinement  levels  desired  (five  levels  for  this  example),  and  b)  the 
minimiim  normal  height  allowed  below  whidi  the  grid  may  not  be  further  refined.  Thus, 
near  the  train  where  the  original  grid  has  fine  resolution,  only  two  or  three  refinement 
levels  were  obtained.  Conversely,  five  refinement  levels  will  be  obtained  at  zones  of  coarse 
initial  grid,  such  as  far  above  the  train. 
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Figure  1-1  shows  the  mesh,  mesh  refinement  levels,  and  pressure  contours  for  the 
complete  computational  domain  and  for  an  expanded  view  around  the  train  at  t=0.  A 
highly  refined  mesh  was  obtained  only  in  the  immediate  vicinity  of  the  shock.  These 
figures  indicate  the  futiUty  of  any  attempt  to  grid  such  a  computational  domain,  which 
is  several  himdred  meters  in  each  direction,  with  a  uniform  fine  resolution  grid  that  will 
yield  captured  shock  wave  thickness  of  less  than  0.6  cm  (about  the  average  shock  thickness 
obtained  in  this  computation). 

A  regular  Mach  reflection  from  the  ramp  is  observed  at  t=0.4  (Figure  1-2).  The  well 
resolved  slip  line  (shown  in  the  density  data,  Figure  l-2b)  demonstrates  one  of  the  advan¬ 
tages  of  employing  density  rather  than  pressure  as  the  critical  refinement  parameter.  A 
comparison  of  the  predicted  pressure  values  behind  the  Mach  stem  with  available  experi¬ 
mental  data  showed  very  good  agreement.  The  incident  shock  hits  the  train  at  t=0.46. 

The  solution  around  the  complete  train  at  t=0.6  (Figure  1-3)  shows  a  curved  shock 
reflection  from  the  top-upstream  comer  of  the  train,  followed  closely  by  an  expansion  fan. 
An  expanded  view  of  pressure  contours  near  the  bottom  of  the  train  at  t=0.5  (Figure  1- 
4a)  shows  the  non-perturbed  incident  shock  (shock  la),  the  main  shock  reflection  from 
the  upstream  surface  of  the  train  (shock  1),  the  rarefaction  waves  expanded  from  the 
upstream  bottom  comer  of  the  train  (wave  5)  and  the  upstream  top  comer  of  the  elevation 
(wave  4)  after  the  main  shock  passage,  the  reflected  ramp  shock  (shock  2),  and  the  curved 
Mach  stem  (shock  3)  that  was  further  weakened  by  the  interaction  with  the  rarefaction 
wave  4,  The  expanded  view  imder  the  train  at  t=0.6  (Figure  l-4b),  shows  the  primary 
reflected  shock  from  the  train  (shock  1),  its  reflection  from  the  top  of  the  elevation  (shock 
la);  the  rarefaction  wave  centered  at  the  upstream  bottom  comer  of  the  train  (wave  5); 
the  reflected  ramp  (shock  2),  its  reflection  from  the  train  (shock  2a),  and  the  rarefaction 
wave  produced  upon  the  impact  of  this  wave  on  the  train’s  upstream  bottom  comer  (wave 
2b);  the  incident  Mach  stem  under  the  train  (shock  3);  and  the  expansion  wave  (wave  4) 
centered  at  the  upstream  comer  of  the  ramp. 

A  comparison  of  the  adapted  mesh  and  mesh  refinement  levels  at  t=0.8  (Figures  1- 
5a  and  l-5b)  with  the  pressure  results  (Figure  l-5c)  demonstrates  the  ability  of  the  grid 
adaptation  routine  to  adapt  to  all  density  gradients.  Due  to  the  initially  finer  grid  under 
the  train,  only  three  refinement  levels  were  obtained  there.  In  contrast,  five  refinement 
levels  were  observed  away  from  the  train,  where  the  initial  grid  was  fairly  coarse.  The 
pressure  results  at  this  time  indicate  that  the  shock  diffraction  over  the  top  upstream 
curved  comer  of  the  train  produced  a  triple  point  connecting  the  incident  shock,  a  Mach 
stem,  and  a  curved  reflected  shock  (Figure  l-5c).  Ramp  shock  reflection  (shock  2  in  Figure 
l-4b)  from  the  front  surface  of  the  train  formed  a  high-amplitude  Mach  stem  traversing 
the  front  siuface  of  the  train. 
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Complex  flow  processes  were  observed  under  the  train.  Examination  of  the  expanded 
views  imder  the  train  at  t=0.7  (Figure  l-4c)  and  t=0.8  (Figure  l-4d)  show  mwy  interesting 
shock  diffraction  processes.  At  t=0.7,  the  incident  Mach  stem  has  just  cleared  the  upstream 
support  beam  (shock  1).  The  reflection  of  the  incident  shock  from  the  beam  has  propagated 
upstream  (shock  3  at  t=0.7  and  shock  la  at  t=0.8)  and,  due  to  its  curvature,  partially 
reflected  from  the  bottom  of  the  train  (shock  lb  at  t=0.8).  The  diffracted  incident  Mach 
stem  (shock  1)  downstream  of  the  beam  has  impacted  on  the  top  of  the  elevation  at  t=0.7, 
and  was  reflected  as  shock  3  at  t=0.8.  Immediately  behind  it,  and  partially  merged,  is 
shock  2d,  which  is  the  transmitted  part  (across  the  beam)  of  shock  2a  at  t=0.6.  The  ramp 
shock  and  its  reflection  from  the  train  (shocks  2  and  2a  in  Figure  l-4b)  are  shocks  2  and 
2a;  the  expanding  rarefaction  waves  2b  and  5  at  t=0.6  (Figure  l-4b)  are  waves  2b  and  5  at 
t=0.7,  and  waves  2c  and  2b,  respectively,  at  t=0.8;  the  reflection  of  the  main  shock  (shock 
1  in  Figure  l-4b)  from  the  top  of  the  elevation  is  shock  la  at  t=0.7,  and  7a  at  t=0.8,  while 
the  reflection  of  shock  2a  at  t=0.6  (which  is  the  reflection  of  the  reflected  ramp  shock  from 
the  train)  is  shock  2c  at  t=0.7,  and  shock  7b  at  t=0.8.  The  evolution  of  these  shocks 
is  important  in  determining  the  vibrations  on  the  train.  Shock  la  at  t=0.7  has  partially 
reflected  from  the  upstream  beam  (shock  5  at  t=0.8,  which  has  also  reflected  from  the 
bottom  of  the  train  as  shock  5a),  and  partially  transmitted  as  shock  4  at  t=:0.8,  which  has 
reflected  from  both  the  bottom  of  the  train  (shock  4a)  and  the  top  of  the  elevation  (shock 
4b).  Similarly,  shock  2c  has  expanded  up  and  downsteam  (shock  6),  and  reflected  from  the 
bottom  of  the  train  (shock  6a).  The  low  pressure  zone  near  the  top  of  the  beam  indicates 
a  high  velocity  zone  created  as  a  result  of  flow  acceleration  aroxmd  the  top  of  the  beam 
and  the  attached  vortex  sheet,  which  has  just  begim  its  roll-up  process. 

While  the  incident  shock  continued  propagating  above  the  train,  the  triple  point  height 
grew  linearly  with  time,  and  the  circular  expansion  waves  from  both  the  bottom  and  top 
upstream  comers  also  continued  growing  linearly  with  time  (Figures  1-3,  l-5c,  l-6a  and  1- 
7c),  the  most  interesting  shock  evolution  processes  occurred  imder  the  train.  Examination 
of  the  pressure  contour  results  tmder  the  train  at  t=0.9  (Figure  l-4e)  and  t=1.0  (Figure 
l-4f)  shows  the  virtual  separation  of  the  shock  system  between  the  support  beams  and  the 
reflected  shocks  upsteam  of  the  beams;  a  separation  resulting  from  the  growing  vortex  near 
the  upstream  support  beam  which  will  accelerate  and  break  any  approaching  shock.  The 
leading  shock  (shock  1)  has  been  completely  merged  with  shock  4  at  t=0.8.  Meanwhile, 
we  observe  the  formation  and  rise  of  the  triple  point  between  shocks  1  and  3,  and  Mach 
stem  Im.  Shock  3  resulted  from  the  merging  of  shocks  3  and  4b  at  t=0.8,  an  almost 
completed  merger  as  two  shocks  are  observed  to  reflect  from  the  bottom  of  the  train  at 
t=0.9  (shocks  3a  and  4b).  Nevertheless,  a  single  reflected  shock  is  observed  after  the  shocks 
were  reflected  again  from  the  top  of  the  elevation  (t=1.0).  Another  shock  system  observed 
between  the  supp>ort  beams  are  shock  6  (at  t=0.9),  which  evolved  from  shock  6a  at  t=0.8. 
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and  its  reflection  from  the  top  of  the  elevation,  shock  6a  (t=1.0) 

Shock  wave  reflection  upstream  is  characterized  both  by  shock  wave  reverberation 
between  the  bottom  of  the  train  and  the  top  of  the  elevation  and  shock  interaction  with 
the  vortex  emanating  from  the  upstream  bottom  comer  of  the  train.  As  an  example 
we  point  to  shock  lb  at  t=0.8  (which  has  just  reflected  from  the  bottom  of  the  train), 
propagating  down,  and  reflecting  from  the  top  of  the  elevation  as  shock  Ic  (t=0.9),  and 
r^ecting  again  from  the  bottom  of  the  train  as  shock  Id  (t=1.0).  Similar  reflection  pattern 
was  observed  for  shock  5a. 

Typical  shock-vortex  interaction  can  be  observed  in  Figures  l-4e  and  l-4f  (t=0.9  and 
1.0,  respectively).  The  reflection  of  the  curved  shock  la  from  the  bottom  of  the  train 
(shock  lb)  resulted  in  the  formation  of  a  triple  point,  with  increased  height  from  the 
surface.  At  t=0.82,  the  Mach  stem  near  the  surface  and  the  two  shocks  (shock  la  and  lb) 
interact  with  the  vortex  attached  to  the  upstream  bottom  comer  of  the  train.  The  result 
is  an  acceleration  and  diffusion  of  the  shock  near  the  surface  of  the  train;  both  effects 
decrease  with  increased  distance  from  the  surface.  These  processes  result  in  the  upstream 
acceleration  and  diffusion  of  compression  wave  Id  between  t=0.8  and  t=1.0.  The  next 
shock  wave  system,  shocks  5  and  5a,  experience  similar  though  weaker  processes,  perhaps 
due  to  the  weakening  of  the  vortex. 

The  very  low  pressure  zone  near  the  upstream  support  beam  indicates  the  rapid  roll¬ 
up  of  a  vortex  sheet  anchored  at  the  top-upstream  comer  of  the  beam  (Figure  l-5c  and 
Figures  l-4e  through  l-4i).  The  vortex  roll-up  process  was  enhanced  by  vortex  interaction 
with  the  system  of  curved  reflected  shocks  reverberating  between  the  bottom  of  the  train 
and  the  top  of  the  elevation.  A  much  weaker  vortex  sheet  has  also  been  formed  near  the 
upstream  bottom  comer  of  the  train.  Mach  number  results  at  t=0.8  and  later  (Figures 
l-6b  and  l-6d  at  t=1.0  and  t=1.2,  respectively)  demonstrate  strong  flow  acceleration  due 
to  the  vortex  roll-up  and  the  reduction  of  the  effective  flow  area  between  the  top  of  the 
elevation  and  the  bottom  of  the  train,  with  the  highest  velocities  above  the  inner-most 
vortex. 

Since  vortex  dynamics  is  a  phenomenon  normally  associated  with  viscous  flow  pro¬ 
cesses,  while  the  present  numerical  scheme  only  solves  the  Euler  equations,  the  source  of 
this  vorticity  deserves  further  discussion. 

Vorticity  production,  as  modeled  in  the  Bjerknes  Theorem  [26],  incorporates  both 
viscous  and  inviscid  contributions,  and  is  expressed  as: 

f  =  -v(i)xv/>+.v^f;  (1.11) 
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where  ^  =  Vorticity;  v  =  Kinematic  Viscosity;  p=  Pressure;  and  p=  Density.  A  time- 
scale  and  order-of-magnitude  analysis  of  these  terms  indicates  that  the  current  vorticity 
production  is  controlled  by  the  baroclinic  effect  (i.e.,  the  first  term  in  the  equation,  which 
represents  the  non-alignment  of  pressure  and  density  gradients),  and  thus  can  be  predicted 
by  the  inviscid  code.  Laser-schilieren  shadowgraph  data  produced  in  a  similar  shock- 
train  interaction  experiment  conducted  using  a  sub-scale  model  in  a  shock  tube  facility 
demonstrated  vortex  sheet  roll-up  within  a  few  microseconds  after  the  shock  diffracted 
around  the  comer,  a  time  period  too  short  for  viscous  processes  to  significsintly  affect  the 
flow.  The  numerical  predictions  for  this  experiment  were  in  excellent  agreement  with  the 
experimental  data.  Due  to  lack  of  space  they  will  be  presented  in  a  future  publication. 

The  incident  shock  propagated  past  the  downstream  support  beam  at  t=1.03.  Re¬ 
sults  at  t=l.l  and  t=1.2  (Figures  l-4g  and  l-4h,  respectively)  show  that  similar  to  the 
shock  diffraction  process  over  the  upstream  beam  (Figure  l-4d  at  t=0.8),  incident  shock 
dif&action  around  the  downstream  beam  resulted  in  the  formation  of  curved  diffracted  and 
reflected  shocks;  shocks  1  and  la  are  the  transmitted  shock  wave  and  its  reflection  from 
the  top  of  the  elevation,  respectively,  while  shocks  2  and  2a  are  the  reflected  shock  and  its 
reflection  from  the  bottom  of  the  train.  The  merged  shock  formed  of  shocks  3  and  6  (at 
t=1.0)  has  been  partially  transmitted  as  shock  6  (t=l.l)  and  was  later  reflected  from  the 
top  of  the  elevation  as  shock  6c  (t=1.2),  and  partially  reflected  from  the  beam  as  shock 
fid  (at  t=l.l  and  t=1.2).  Shocks  4  and  6a  at  t=1.0  continued  to  reverberate  and  produce 
shocks  4,  6a  and  fib  at  t=l.l,  and  shocks  4,  4a  and  4b  (beam  reflection),  and  fib  at  t=1.2. 
The  vorticity  contour  plots.  Figures  l-6c  and  l-6e,  at  t=1.0  and  t=1.2,  respectively,  indi¬ 
cate  significant  enhancement  of  the  rolled-up  vortex  sheet  behind  the  upstream  beam  due 
to  the  interaction  with  the  family  of  curved  reflected  shocks,  and  the  initiation  of  a  vortex 
sheet  behind  the  downstream  beam. 

The  successful  continued  grid  adaptation  is  demonstrated  again  at  t=1.4.  A  com¬ 
parison  of  the  computational  grid  (Figure  l-7a),  grid  refinement  level  (Figure  l-7b),  and 
pressure,  Mach  niimber  and  vorticity  contours  (Figiires  l-7c,  l-7d  and  l-7e,  resjjectively) 
demonstrate  the  ability  of  the  new  methodology  to  adapt  to  both  strong  and  weak  density 
gradients. 

The  shocks  that  propagated  imder  the  train  emerged  at  t=1.22.  The  higher  stagnation 
pressure  tmder  the  train  resulted  in  faster  propagation  speed  than  above,  and  thus  the 
incident  shock  propagating  below  has  emerged  a  short  time  earlier.  Pressure  results  at 
t=1.4  (Figure  l-4i)  indicate  that  four  primary  shocks  emerged:  the  incident  shock  (shock 
1);  its  reflection  from  the  top  of  the  elevation  (shock  la);  shock  6c;  and  shock  6,  which  has 
been  weakened  by  its  upward  expansion  and  broken  due  to  its  interaction  with  the  stronger 
shocks  la  and  6c.  The  expansion  of  these  shocks  aroimd  the  bottom  downstream  comer  of 


11 


the  train  and  the  comer  of  the  elevation  resulted  in  the  immediate  formation  of  a  counter¬ 
clockwise  rotating  vortex  sheet  anchored  at  the  comer  and  a  cylindrical  expansion  wave. 
The  bottom-produced  expansion  could  not  prop>agate  upstream  beyond  the  downstream 
beam,  due  to  the  transonic  flow  between  the  beam  and  the  bottom  of  the  train. 

Among  the  upstream  reflected  shocks  observed  are  the  primary  reflection  from  the 
train  upstream  surface,  reflections  of  the  Mach  stem  from  the  train,  and  several  shock 
reflections  from  the  upstream  support  beam.  The  system  of  shocks  between  the  beams 
continued  reverberating  and  producing  new  shocks.  The  interactions  between  these  curved 
shocks  with  the  rolled-up  vortex  attached  to  the  upstream  beam  resulted  in  the  break-up 
of  this  sheet,  as  observed  at  t=1.4  (Figure  l-7e)  and  t=1.6  (Figure  l-8a). 

Figure  1-8  shows  vorticity  contours  at  several  times.  While  vorticity  production  was 
fairly  low  for  shock  diffraction  over  the  rounded  upper  downstream  comer  of  the  train, 
the  system  of  shock  waves  reverberating  between  the  support  beams  and  between  the  top 
of  the  elevation  and  the  bottom  of  the  train  both  enhanced  and  broke  up  the  rolled-up 
vortices  attached  to  the  support  beams,  as  shown  in  Figures  l-8a,  l-8b,  l-8c,  l-8d  and  l-8e 
at  t=1.6, 1.8,  2.0, 2.2  and  2.4,  respectively.  Meanwhile,  the  vortices  attached  to  the  bottom 
comers  of  the  train  continued  their  roll-up  process;  the  one  attached  to  the  downstream 
comer  rolled  up  about  four  times  at  t=2.4  (Figtire  l-8e). 

The  shocks  diffracting  arotmd  the  downstream  side  of  the  train  interacted  with  the 
opposite  comers  at  t=2.6.  The  computational  grid,  the  grid  adaptation  levels,  and  the 
pressure  and  Mach  number  contour  resiilts  at  t=2.4  (Figur«;s  l-9a  through  l-9d)  demon¬ 
strate  the  excellent  shock  adaptation  capabilities  of  the  developed  methodology,  even  after 
propagating  long  distances.  The  large  disparity  between  the  long  shock  propagation  dis¬ 
tance  and  the  small  size  of  the  resolved  vortices  or  the  captured  shocks  thickness  clearly 
demonstrates  the  economical  advantage  of  applying  an  adaptive  refinement /coarsening 
scheme  as  compared  to  a  fixed-mesh  scheme. 

The  shock  that  has  diffracted  over  the  top  interacted  at  t=2.6  with  the  vortex  sheet 
anchored  at  the  bottom  downstream  comer  (Figure  l-8f).  Since  this  curved  diffracted 
shock  carried  vorticity  of  opposite  sign  and  of  about  equal  magnitude  to  the  vorticity  of  the 
attached  vortex  sheet,  the  interaction  caused  the  eventutil  detachment  of  the  vortex  sheet 
from  the  comer  and  its  downstream  convection.  Flow  area  constriction  due  to  the  vortices 
shed  from  both  the  upstream  and  downstream  beams  (Figures  l-7d  and  l-9d)  resulted 
in  flow  acceleration  to  sup>ersonic  velocities,  specifically,  downstream  of  the  downstream 
beam  (M=1.28).  The  instantaneous  local  flow  area  expansion  near  the  downstream  bottom 
comer  of  the  train  resulted  in  the  formation  of  a  normal  shock  attached  to  that  comer,  and 
a  deceleration  to  M=0.65.  Further  downstream,  effective  flow-aiea  constriction  between  the 
vortex  attached  to  this  comer  and  the  vortex  attached  to  the  downstream  beam  accelerated 


the  flow  to  M=1.2,  which  was  further  followed  by  an  oblique  shock  as  the  effective  flow-area 
increased  again. 

The  downward-propagating  shock  downstream  of  the  train  hit  the  bottom  bevel  at 
t=3.1  and  reflected  upward.  Meanwhile,  the  shock  that  emerged  from  under  the  train 
continued  its  downstream-and-up  propagation,  though  shock  propagation  upward  was  in¬ 
hibited  by  the  transonic  flow  above  the  train.  This  blockage  effect  temporarily  trapped  a 
relatively  high  pressure  zone  downstream  of  the  train,  which  could  not  be  relieved  through 
the  channel  imder  the  train  due  to  the  transonic  flow  in  this  zone.  Meanwhile,  the  vor¬ 
tices  shed  from  the  downstream  beam  and  bottom  comer  were  convected  downstream  by 
the  flow,  as  shown  in  Figures  l-8g,  l-8h,  l-8i,  l-8j  and  l-8k,  at  t=3.2,  3.6,  4.2,  5.0  and 
6.0,  respectively.  These  results  also  show  the  break-up  of  the  vortex  sheet  attached  to 
downstream  beam,  the  roll-up  of  the  broken  sheet  into  a  weedc,  clockwise-rotating  vortex 
(Figxnre  l-8h),  and  the  rolling  of  the  pair  of  counter-rotating  vortices  by  the  flow  (Figures 
l-8i,  l-8j  and  l-8k). 

The  computation  was  continued  to  t=8.0,  long  after  the  system  of  shocks  cleared  the 
train.  At  this  time  the  flow  above  the  train  was  still  transonic  after  reaching  a  maximum 
Mach  number  of  M  =  1.35  at  t=4.5.  Thus,  the  upward  propagation  of  the  shock  reflected 
from  the  downstream  bevel  was  delayed,  and  its  influence  wa'^  still  felt  on  the  top  of  the 
train  at  t=6.0  (Figure  1-10).  All  shocks  had  propagated  very  long  distances,  but  were 
still  captured  as  sharp,  nonoscillatoiy  discontinuities.  The  captured  shocks’  thickness  has 
not  changed  over  the  past  several  thousand  steps,  as  all  shocks  were  still  captured  over 
two  to  three  elements,  while  the  element  size  aroimd  the  shock  has  not  changed  due  to 
the  adaptive  procedure.  These  results  demonstrate,  at  least  for  the  present  application, 
the  advantage  of  adaptive-mesh  schemes  over  fixed-grid  schemes.  Achieving  similar  shock 
resolution  with  a  fixed-mesh  methodology  would  have  required  maintaining  a  fine  grid  vir¬ 
tually  everywhere  in  the  computational  domain,  at  a  significantly  (perhaps  prohibitively) 
higher  computational  cost.  Grid  refinement  was  observed  as  the  shocks  propagated  into 
new  areas,  while  grid  coarsening  was  observed  in  areas  already  traversed  by  the  shocks. 
Finally,  it  should  be  noted  that  the  solution  obtsdned  above  the  train  at  later  times  in¬ 
dicates  a  compression  wave  rather  than  a  shock.  This  is  the  correct  physical  solution  as 
the  initial  reflected  shock  above  the  train  (Figure  l-5a  at  t=l,0  and  Figure  l-7a  at  t=1.4) 
was  significantly  weakened  by  the  strong  expansion  over  the  train  (Figure  l-9c  at  t=2.4). 
Results  obtained  in  other  simulations  using  significantly  stronger  shocks  demonstrated  a 
reflected  shock  above  the  trmn,  even  at  very  late  times. 

Comparisons  of  experimental  pressure  and  impulse  time  histories  and  the  correspond¬ 
ing  numerical  predictions  are  shown  in  Figures  1-11  at  selected  stations  aroimd  the  train. 
Stations  a  and  b  were  located  on  the  upstream  surface  of  the  train,  station  c  on  top. 


13 


stations  d  and  e  on  the  downstream  surface  and  stations  f  and  g  on  the  bottom.  The 
experimental  data  were  obtained  in  a  large  shock  tube  facility.  Since  the  experimental 
model  spanned  a  significant  portion  of  the  width  of  the  shock  tube,  it  was  appropriate  to 
model  this  flow  with  a  two-dimensional  algorithm.  The  numerical  pressure-time  histories 
and  impulses  (temporally  integrated  pressures)  are  shown  as  thicker  lines  to  t=8.0.  The 
experimental  pressure  and  impulse  data  are  shown  in  thinner  lines  to  t=10.0.  Very  good 
agreement  is  demonstrated  at  all  locations.  Physical  processes  controlled  by  shock  physics 
were  predicted  extremely  well.  For  instance,  very  good  agreement  is  demonstrated  between 
the  predicted  and  measured  presstire-time  histories  at  stations  six  and  seven,  located  un¬ 
der  the  train;  these  stations  exhibit  a  multi-shock  system  that  resulted  from  the  repetitive 
shock  reflection  from  the  top  of  the  elevation,  the  support  beams,  and  the  bottom  of  the 
train.  In  contrast,  the  late-time  experimental  data  (well  past  the  diffraction  phase)  for  sta¬ 
tions  on  the  train  top  and  back  indicate  periodic  3-D  vortex  shedding,  a  physical  process 
that  cannot  be  properly  modeled  by  the  present  2-D  model. 

A  final  note  relating  to  the  computationsJ  performance  of  the  code:  computation 
time  for  approximately  6800  time  steps  consumed  approximately  six  hours  of  CPU  time 
on  a  Cray  2/8-128  computer  (single  processor).  Computation  time  was  approximately  60 
microseconds  per  node  per  time  step. 

1.4  SUMMARY  AND  CONCLUSIONS. 

A  new  transient,  two-dimensional,  finite-element  shock  capttiring  scheme  on  unstruc¬ 
tured  grids  was  applied  to  the  study  of  shock  interaction  with  a  train  suspended  above  a 
rigid  elevated  surface.  The  area  between  the  train  and  the  surface  was  partially  blocked 
by  the  train  support  beams,  resulting  in  complex  shock  diffraction  processes.  The  results 
demonstrate  the  ability  of  the  new  adaptive  refinement /coarsening  algorithm  to  resolve 
shocks  in  a  sharp-nonoscillatory  manner.  In  addition  to  interesting  shock  wave  propaga¬ 
tion  and  interaction  processes,  the  results  demonstrate  the  capability  of  the  new  code  to 
capture,  and  define  in  great  detail,  vortices  shed  from  sharp  comers. 

Among  the  more  interesting  shock  propagation  processes  observed  were  shock  diffrac¬ 
tion  aroimd  sharp  comers  and  the  immediate  formation  of  vortices,  no  doubt  due  to  the 
baroclinic  effect;  interaction  of  diffracted  curved  shocks  with  vortex  sheets  carrying  vor- 
ticity  of  identical  or  opposite  sign,  resulting  in  either  vortex  enhancement  or  break-up; 
containment  of  the  high  pressure  zone  (i.e.,  blockage  of  upward  propagation  of  shocks)  on 
the  downstream  side  of  the  train  by  the  supersonic/transonic  flow  above  the  train;  and 
the  blockage  of  the  upstream  propagation  of  the  rarefaction  wave  in  the  channel  imder  the 
trmn  due  to  the  transonic  flow  above  the  downstream  support  bezun. 
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SECTION  2 

NUMERICAL  DESIGN  OF  A  PASSIVE  SHOCK  REFLECTOR 


A  recently  developed  transient,  two-din^sional  finite-element  shock-capturing  scheme 
on  unstructured  grids,  was  applied  to  the  study  of  a  shock  interacting  with  a  passive  shock 
reflector  placed  at  an  opening  of  a  vented  enclostire.  The  objective  of  this  effort  was  to 
design  a  passive  device  that,  while  allowing  the  ventilation  of  the  enclosure  under  steady 
conditions,  will  prevent  blast-waves  impinging  on  the  wall  firom  entering  the  enclosure  when 
the  structure  is  impacted  by  a  shock.  A  passive  reflector  (‘*grill”)  was  designed  to  take 
advantage  of  two  basic  fluid  dynamic  processes:  shock  dif&action  around  sharp  comers, 
and  vortex  shedding  firom  sharp  comers.  As  the  shock  propagated  through  the  grill,  shock 
diffraction  around  several  comers  resulted  in  an  order-cf-magnitude  shock  amplitude  re¬ 
duction,  while  transforming  the  steep-fionted  shock  into  a  compression  wave.  The  vortex 
shedding  process  significantly  reduced  the  effective  flow  area,  and  limited  the  quasi-steady 
flow  behind  the  shock  from  entering  the  enclosure.  F^om  the  numerical  point  of  view, 
the  results  demonstrate  the  capability  of  the  developed  adaptive  refinement /coarsening 
algorithm  to  properly  adapt  to  weak  shocks,  rarefaction  waves  and  other  weak  flow  gra¬ 
dients,  and  the  resultant  excellent  resolution  of  the  captured  flow  features.  In  addition 
to  interesting  shock  diffraction  and  propagation  phenomena,  the  results  demonstrate  the 
capability  of  the  new  code  to  capture,  and  define  in  great  detail,  vortex  sheets  shed  from 
sharp  comers.  Finally,  the  results  demonstrate  that  when  a  design  project  is  based  on  the 
imderstanding  and  application  of  basic  fluid-dynamic  mechanisms,  it  is  possible  to  obtain 
significant  improvements  in  performance,  specifically,  an  order-of-magnitude  reduction  of 
shock  amplitude  entering  the  enclosure. 

2.1  INTRODUCTION. 

During  the  past  ten  years,  the  CFD  community  has  experienced  a  proliferation  of  a 
multitude  of  shock  captviring  schemes  whose  ultimate  objective  is  the  sharp,  nonoscillatory 
capttuing  of  transient  shocks,  even  after  propagating  long  distances  [1-3].  The  practical 
objective  of  these  methodologies  is  to  predict  the  loads  exerted  by  shocks,  initiated  by 
nuclear  or  High-Explosive  (HE)  bursts,  on  stationary  or  moving  stmctures  located  a  long 
distance  from  the  burst  point.  The  results  of  such  simulations  should  help  improve  target¬ 
ing  specification  requirements  (for  offensive  purposes)  or  hardening  criteria  (for  defensive 
objectives).  FVom  the  numerical  algorithm  development  point  of  view,  the  large  propa¬ 
gation  times  and  distances  that  the  shocks  are  likely  to  traverse  fiom  burst  to  impact 
point  pose  significantly  greater  demands  than  required  from  traditional  shock  capttuing 
schemes.  These  schemes  are  only  intended  to  produce  steady,  sharp,  nonoscillatory  shocks 
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on  airfoils,  etc.  [7]  to  yield  the  steady  lift  and  drag  forces.  Past  simulations  of  shock  wave 
propagation  over  long  distances  relied  on  either  the  fixed-mesh  structured  grid  approach 
[1,2]  or  sliding  refined  zones  [20],  intended  to  continuously  surround  the  shocks  with  a  finer 
mesh  than  elsewhere  in  the  domain.  The  first  approach  yielded  coarser,  diffused  shocks 
after  some  propagation  distance,  since  even  for  two-dimensional  calculations  it  is  not  eco¬ 
nomically  feasible,  even  with  present  day  class  VI  computers,  to  grid  a  domain  of  several 
himdreds  or  thousands  of  meters  with  1  cm  size  cell.  Worse  yet,  in  this  investigation  it 
is  necessary  to  compute  flow  diffraction  around  a  grill  with  small  gaps,  which  requires 
mUlimeter-scaU  resolution.  The  second  approach  was  much  more  economical,  as  only  part 
of  the  domain  was  finely  refined.  Nevertheless,  the  computational  resources  required  for 
two-dimensional  computations  (both  memory  and  computational  time)  were  very  high  and 
required  knowledge  of  the  primary  shock  location,  while  three-dimensional  computations 
would  have  been  prohibitively  expensive  for  production-type  runs.  Therefore,  the  approach 
employed  here  was  a  natural  evolution  of  past  deficiencies,  i.e.,  employ  a  numerical  scheme 
capable  of  d3rnamically  adapting  to  traveling  shocks  and  other  fiow  discontinuities.  This 
approach  allows  us  to  use  coarse  grid  resolution  everywhere  in  the  computational  domain 
where  fiow  gradients  are  low,  while  appl}ring  very  fine  grid  resolution  wherever  better  flow 
resolution  is  required.  Hence,  shocks  should  always  advance  into  finely  refined  zones,  while 
grid  coarsening  should  be  obtained  in  areas  already  traversed  by  the  shocks. 

The  objective  of  this  study  was  to  use  the  recently  developed  adaptive  grid  method¬ 
ology  on  unstructured  grids  [11]  to  evaluate  the  effectiveness  of  candidate  passive  shock 
reflector  devices.  This  device  was  placed  at  an  opening  in  the  wall  of  a  semi-enclosed 
room.  A  radiator,  or  another  heat-exchange  device,  was  placed  adjacent  to  the  grill  within 
the  room.  The  shock  reflector  is  supposed  to  allow  free  exchange  of  air  during  normal 
operations,  while  preventing  blast  waves  from  entering  the  enclosure  when  the  structure 
is  impacted  by  a  shock.  A  simplified  schematic  of  the  structure  and  the  shock  reflector 
(often  referred  to  as  “grill” )  com|>osed  of  simple,  readily  available,  90“  chevrons  is  shown 
in  Figure  2-1.  Three  layers  of  chevrons  were  employed  in  this  specific  design.  A  series  of 
flat-plates  was  placed  in  front  of  the  chevrons  to  block  the  flow  from  entering  the  spacing 
between  the  chevrons. 

2.2  NUMERICAL  RESULTS. 

The  transient,  two-dimensional  finite-element  shock-capturing  adaptive  scheme  on 
unstructured  grids  described  in  Section  1  (FEFL027)  was  applied  to  the  study  of  a  shock 
interaction  with  a  passive  shock  reflector  placed  at  an  opening  of  a  vented  enclosure  (Figure 
2-1).  The  computational  parameters  used  in  the  simulation  included: 
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•  Number  of  time  steps  per  refinement:  seven 

•  Niimber  of  refinement  levels:  five 

•  Minimum  normal  size  allowed:  25  /xm 

•  Key  variable  for  the  error  indicator:  density 

•  Courant  number:  0.4 

The  color  contour  results  have  256  contovirs,  where  blue  represents  the  lowest  instantaneous 
local  value  and  magenta  the  highest.  The  color  designations  for  the  mesh  refinement  levels 
are  blue  for  the  original  mesh,  and  green,  yellow,  red,  magenta  and  cyan  for  one,  two,  three, 
four  and  five  refinement  levels,  respectively:  one  level  of  mesh  refinement  is  defined  as  the 
division  of  a  trian^e  into  four,  performed  by  connecting  the  mid-faces.  Two  parameters 
controlled  mesh  refinement:  a)  the  maximum  number  of  levels  desired:  five  levels  for  this 
simulation;  and  b)  the  minimum  normal  height  allowed:  here,  25  pm,  compared  with 
the  length  of  the  computational  domain  (approximately  five  meters).  A  similarly  resolved 
fixed-mesh  methodology  would  require  approximately  1.2*  10**’  nodes  -  an  impossible  task. 

The  initial  conditions  and  the  geometric  details  of  this  problem  are  shown  in  Figure 
2-1,  which  shows  both  the  initial  pressure  contours  for  the  complete  computational  domain 
and  for  an  expanded  view  near  the  middle  of  the  griU.  The  incident  weak  shock,  a  step 
function  with  an  amplitude  of  10  psi  overpressure  (a  value  chosen  for  demonstration  pur¬ 
poses  only),  was  placed  at  tsQ  to  the  left  of  the  enclosed  chamber.  The  shock  propagated 
from  left  to  right.  The  flow  behind  the  shock  was  locally  subsonic,  while  the  shock  Madi 
ntimber  was  about  1.4.  Figure  2-1  actually  shows  2-D  side-view  schematics  of  the  structure 
and  the  shock  reflector  (which  will  be  referred  to  henceforth  as  the  ‘*grill”).  Three  layers 
of  simple,  readily  available  90**  chevrons  were  employed  in  this  specific  design.  An  array 
of  bars  (flat  plates)  was  placed  in  front  of  the  chevrons  to  block  the  flow  from  entering 
the  spacing  between  the  chevrons.  This  design  was  intended  to:  a)  force  the  shock  pene¬ 
trating  between  the  chevrons  to  diffract  at  large  angles  to  optimize  shock  wave  damping 
by  the  rarefaction  waves  originating  at  the  comers;  b)  optimize  shock  deflection  on  direct 
impact  at  any  angle;  and  c)  produce  vortidty  behind  the  sharp  comers  that  will  reduce 
the  effective  flow  area  through  the  chddng  of  the  flow  between  the  chevrons. 

In  the  following  analysis  of  the  results,  we  will  initially  describe  the  large-scale,  mul¬ 
tiple  shock-diffnu:tion  processes  from  initiation  to  the  end  of  the  simulation;  we  will  then 
re-examine  flow  expansion  aroimd  the  top  of  the  bar,  and  conclude  with  the  discussion  of 
the  flow  aroimd  the  upstream  comer  of  the  first  chevron.  To  enhance  understanding  of 
the  physical  processes  and  better  exhibit  the  micro-scale  physics,  rather  than  show  flow 
results  for  the  complete  computational  domain,  we  will  show  expanded  views  of  only  one 
row  of  bar  and  chevrons  combination  located  near  the  center  of  the  grill;  the  lower  quarter 
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of  the  grill  shown  in  Figure  2- lb.  Thus,  both  the  top  and  the  bottom  boundaries  of  this 
domain  are  actually  planes  of  symmetry.  For  convenience  sake,  they  will  be  referred  to  in 
the  following  discussion  as  the  top  and  bottom  planes,  while  waves  transmitted  from  the 
other  side  of  the  boundary  will  be  referred  to  as  reflected  waves. 

Figures  2-1  through  2-11  show  combinations  of  computational  mesh,  grid  refinement 
levels,  pressure,  vorticity  and  Mach  number  contours  at  several  times.  Figure  2-2  shows 
the  mesh  refinement  leveb,  density  and  Mach  number  contoiurs  at  t=0.125  ms.  To  better 
appreciate  the  dimensions  of  the  physical  processes  simulated  and  the  accuracy  of  the  nu¬ 
merical  scheme,  it  is  noted  that  the  thickness  of  the  bar  was  1  cm,  or  10,000  ftm.  Excellent 
mesh  adaptation  to  the  incident  and  reflected  shocks,  as  well  as  the  large  recirculation 
areas  on  top  of  the  bar,  are  demonstrated.  The  Mach  number  of  the  flow  accelerating 
around  the  comers  is  0.95,  increasing  with  time.  At  t=0.15  ms,  the  flow  Mach  number 
above  the  bar  was  already  1.05.  The  importance  of  this  fact  will  be  highlighted  later  in 
the  discussion. 

Rarefaction  wave  expansion  from  both  comers  is  clearly  shown  in  both  Figures  2- 
2a-c  (t=0.125  ms),  and  Figure  2-2d  (t=0.18  ms),  at  which  time  the  rarefaction  wave 
originated  at  the  upstream  comer  of  the  bar  reflected  from  the  bottom  plane  (or,  since 
this  is  really  a  plane  of  symmetry,  the  arrival  of  the  rarefaction  wave  from  the  opposite  side 
of  the  bar).  This  process  is  vividly  demonstrated  by  the  transition  of  the  planar  reflection 
(Figure  2-2c)  to  a  curved  reflected  shock  (Figure  2-3b).  The  incident  shock  impinged 
on  the  first  chevron  at  ts:0.21  ms.  Simultaneoiisly,  the  curved,  incident  shock  that  has 
diffracted  arotmd  the  downstream  comer  of  the  bar,  has  reflected  from  the  bottom  plane. 
The  upstream-propagated  initial  reflection  from  the  bar  has  reflected  from  the  top  plane 
(actusdly,  the  reflected  wave  from  the  adjacent  bar  above).  The  Mach  number  resiilts 
(Figure  2-3b)  show  that  the  flow  above  the  bar  has  reached  M=1.25.  These  results  also 
show  the  merging  of  the  two  recirculating  zones  above  the  bzu-,  incident  shock  reflection 
from  the  upstream  front  of  the  first  chevron,  the  instantaneous  initiation  of  rarefaction 
waves  at  both  comers,  and  the  simple  shock  reflection  produced  as  the  incident  shock 
propagated  toward  the  convex  center  of  the  first  chevron  (an  obvious  wave  focusing  point). 
Figure  2-3c  (t=:0.245  ms)  shows  the  continued  evolution  of  these  processes:  the  incident 
shock  has  diffracted  around  the  upstream  comer  of  the  first  chevron;  two  rarefaction  waves 
have  emanated  from  the  two  upstream  comers  of  the  first  chevron,  the  one  from  the  bottom 
comer  was  fairly  cylindrical,  while  the  one  from  the  top  comer  was  contained  within  the 
reflected  flattened  shock.  Incident  shock  wave  diffraction  over  the  90°  comers  of  the  front 
bar,  followed  by  its  diffraction  around  the  135°  comer  of  the  first  chevron,  resulted  in 
significant  pressture  amplitude  reduction  of  the  shock  penetrating  past  the  first  chevron  - 
only  about  3.5  psi  at  this  time. 
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The  density  and  Mach  number  restilts  at  t=0.27  ms,  shown  in  Figures  2-3d  and  e, 
respectively,  show  the  continued  growth  of  the  merged  vortices  at  the  top  of  the  front  bar, 
where  the  single  vortex  was  anchored  at  the  upstream  comer;  the  reflection  of  the  shock 
that  originally  reflected  from  the  first  chevron,  from  the  bottom  plane;  and  the  eminent 
interaction  of  the  incident-shock  reflection  from  the  bottom  plane  with  the  two  vortices 
(the  one  anchored  at  the  top  upstream  comer  of  the  bar  and  the  one  anchored  at  the 
upstream  bottom  comer  of  the  first  chevron). 

Shock  wave  focusing  at  the  convex  center  of  the  first  chevron  occurred  at  t=0.315  ms. 
The  computational  grid,  mesh  refinement  levels,  density  and  Maurh  number  contours  at 
t=0.35  ms,  shown  in  Figures  2-4a-d,  respectively,  demonstrate  the  following:  the  excellent 
shock  adaptation;  the  very  small  ratio  of  area  well  refined  (i.e.,  occupied  by  the  shocks) 
to  the  complete  computational  domain;  the  supersonic  flow  over  the  first  layer  of  chevrons 
(M=1.24  at  this  time);  the  successful  mesh  coarsening  after  shock  passage;  and  finally, 
the  large  amplitude  of  the  reflected  shock  from  the  center  of  the  first  chevron  (about  65 
psi),  much  higher  than  would  have  been  obtained  for  a  regular  shock  reflection.  The  large 
amplitude  reflected  shock  will  be  shown  at  later  times  (Figures  2-5  to  2-8)  to  create  a  high- 
pressure  reflected  jet  that  will  persist  for  a  very  long  time.  Incidently,  from  a  practical 
point  of  view,  the  large  pressure  will  not  damage  this  chevron  as  the  pressure  decays 
exponentially  within  2-3  ms,  and  hence,  the  integrated  impulse  is  below  bending  failure 
level.  Both  the  mesh  refinement  and  density  contours  show  systems  of  weak-amplitude 
shocks  (shocklets)  emanating  from  the  large  separation  zones  near  the  comers.  We  will 
expand  on  the  importance  of  this  subject  later  in  the  discussion. 

The  downstream-propagated  weakened  incident  shock  expanded  over  the  second  layers 
of  chevrons  at  t=0.45  ms  (Figure  2-5a).  While  the  incident  shocks  focused  at  the  convex 
center  of  the  second  chevron  at  t=0.42  ms,  forming  another  fairly-high  amplitude  reflected 
wave  (about  18  psi  at  t=0.45  ms)  that  will  help  in  blocking  the  incoming  flow,  the  incident 
shock  expanded  over  another  pair  of  90*  comers,  and  was  further  weakened  by  another 
pair  of  cylindrical  rarefaction  waves  expanded  from  both  comers  of  the  second  chevron 
(especially  the  downstream  comer).  Generally,  the  diffraction  processes  around  the  second 
and  third  layers  of  chevrons  were  very  similar  to  the  diffraction  around  the  first  chevron, 
i.e.,  shock  wave  focusing  toward  the  center  of  the  chevron,  instantaneous  localized  pressure 
increase,  and  a  much  l2U'ger  amplitude  (due  to  wave  focusing)  reflected  shock.  Nevertheless, 
the  incident  and  reflected  pressure  amplitudes  decreased  drastically  from  layer  to  layer 
(about  65  psi  for  the  first,  22  psi  for  the  second  zmd  about  4.5  psi  for  the  third),  numbers 

The  incident  wave  impacted  on  the  third  layer  of  chevrons  at  about  t=0.47  ms,  and  on 
its  middle  at  t=0.52  ms.  At  t=0.50  ms  (Figures  2-5b  and  c),  the  amplitude  of  the  incident 
shock  about  to  impact  on  the  third  chevron  was  only  about  2.2  psi,  resulting  from  the 
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iniiltiple  diffraction  processes  it  had  experienced.  The  shock  that  has  dif&acted  around 
the  third  layer  of  chevrons  was  further  weakened  by  another  135°  diffraction  process.  In 
fact,  energy  dissipation  of  the  incident  shock  by  interaction  with  eight  different  rarefactions 
waves  (i.e.,  propagating  past  eight  comers),  resulted  in  very  strong  attenuation  of  the 
incident  10  psi  shock  wave  to  a  mere  compression  wave,  with  sin  amplitude  of  about  1.7 
psi.  The  results  at  this  time  (Figure  2-5c)  also  show  the  system  of  shocks  emanating  from 
the  complex  separated  vortex  downstream  of  the  front  bar.  This  multiple-shock  system 
resulted  from  the  interaction  between  the  mean  supersonic  flow  between  the  bar  and  the 
flrst  chevrons,  and  the  convoluting  vortex.  This  shock  system  evolves  in  time  as  the  vortex 
breaks  down. 

The  vortex  roll-up  process  was  enhanced  by  vortex  interaction  with  the  system  of 
curved  reflected  shocks  reverberating  between  the  bar  and  the  first  chevron,  and  similarly, 
at  other  times,  between  the  various  chevrons.  Since  vortex  dynamics  is  a  phenomenon 
normally  associated  with  viscous  flow  processes  while  the  present  numerical  scheme  only 
solves  the  Euler  equations,  the  sovirce  of  this  vorticity  must  be  addressed. 

Vorticity  production,  as  modeled  in  the  Bjerknes  Theorem  [24],  incorporates  both 
viscous  and  inviscid  contributions,  and  is  expressed  as: 

|[  =  -v(l)xVP  +  .V^{.  (2.1) 

where  ^  —  Vorticity;  v  =  Kinematic  Viscosity;  P  =  Pressure;  and  p  —  Density.  A  time- 
scale  and  order-of-magnitude  walysis  of  these  terms  indicates  that  the  current  vorticity 
production  is  controlled  by  the  baroclinic  effect  (i.e.,  the  first  term  in  the  equation  that 
represents  the  non-alignment  of  pressure  and  density  gradients),  and  thus  can  be  predicted 
by  the  inviscid  code.  Past  experience  with  shock-train  interaction  [25]  demonstrated  vortex 
sheet  roll-up  within  a  few  microseconds  after  the  shock  diffracted  around  the  comer,  a 
time  period  too  short  for  viscous  processes  to  significantly  affect  the  flow.  It  must  be 
noted  that  we  are  fully  aware  of  the  contribution  of  numerical  viscosity  to  the  vortex 
production  process.  Nevertheless,  based  on  past  experience  with  similar  computations  and 
comparisons  to  experimental  data,  which  were  in  excellent  agreement  [26],  we  believe  this 
effect  to  be  of  secondary  importance. 

The  computational  mesh,  mesh  refinement  levels,  density  and  Mach  number  contours 
at  t=0.6  ms,  are  shown  in  Figures  2-6a-d,  respectively.  These  results,  and  the  density 
results  shown  later  at  t=0.65  ms  (Figure  2-7),  demonstrate  that  the  wave  entering  the 
enclosed  chamber  is  a  weakened  compression  wave,  not  a  steep-  fronted  shock  wave.  The 
distinction  is  important  not  only  because  of  the  mnplitude  difference,  but  because  the  two 
waves  would  stagnate  to  vastly  different  pressure  values.  In  addition,  structural  response 
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to  the  instantaneous  loading  of  a  shock  is  significantly  more  severe  than  the  slower  loading 
of  a  compression  wave.  Hence,  if  a  heat  exchange  device  would  have  been  placed  beyond 
the  third  layer  of  chevrons,  the  difference  between  a  shock  and  a  compression  wave  would 
have  meant  the  difference  between  system  survival  and  destruction.  In  addition,  this  figure 
shows,  in  great  detail,  the  excellent  adaptation  of  the  multiple  shock-wave  system  and  the 
sharpness  of  the  shocks.  The  refinement  level  tolerance  for  this  calculation  was  set  at  a 
very  low  level.  Therefore,  even  the  weak  compression  waves  entering  the  chamber  are  still 
captured  (at  least  one  level  of  refinement). 

The  computation  was  continued  to  t=5.0  ms,  long  after  the  system  of  shocks  cleared 
the  grill.  Only  four  levels  of  refinement  were  used  at  this  time,  with  a  minimum  nor¬ 
mal  of  0.2  mm  (200  pm).  The  results  at  t=1.9  ms  (Figure  2-8)  demonstrate  that:  a)  All 
shocks  had  propagated  very  long  distances,  but  were  still  captured  as  sharp,  nonoscillatory 
discontinuities;  and  b)  The  captured  shocks’  thickness  did  not  change  over  several  thou¬ 
sand  steps,  as  all  shocks  were  still  captured  over  two  to  three  elements,  while  the  element 
size  around  the  shock  remained  vmchanged,  due  to  the  adaptive  procedure.  These  results 
demonstrate,  at  least  for  the  present  application,  the  advantage  of  adaptive-mesh  schemes 
over  fixed-grid  schemes.  Achieving  similar  shock  resolution  with  a  fixed-mesh  methodology 
would  have  required  maintaining  a  fine  grid  virtually  everywhere  in  the  computational  do¬ 
main,  at  a  significantly  (perhaps  prohibitively)  higher  computational  cost.  Grid  refinement 
was  observed  as  the  shocks  propagated  into  new  areas,  while  grid  coarsening  was  observed 
in  areas  already  traversed  by  the  shocks;  a)  The  system  of  compression  waves  entering  the 
enclosure  had  steepened,  as  shoxild  be  expected  from  traveling,  finite-amplitude  compres¬ 
sion  waves;  b)  While  the  pressure  outside  the  chevrons  is  about  12  psi  (red  levels  at  this 
figure),  the  pressiue  inside  was  less  than  1.5  psi;  c)  Vortices  are  observed  to  shed  firom  all 
sharp  comers;  and  d)  The  vortices  shed  from  the  first  and  second  layers  of  chevrons  have 
just  about  blocked  the  entire  area  between  the  first  and  second  chevrons.  Mach  number 
contours  data  demonstrate  that  the  flow  through  the  passage  between  the  chevrons  in 
the  first  layers  had  accelerated  to  supersonic  speeds  and  formed  a  normal  shock  ahead 
of  the  center  of  the  second  layer  chevron.  From  there,  the  flow  has  to  turn  90®,  literally 
squeezing  between  vortices  shed  from  the  first  and  second  layers  of  chevrons,  accelerate 
to  supersonic  speeds  again,  and  then  turn  again  by  about  135®  to  pass  through  the  next 
opening  restricted  by  the  vortex  shed  from  the  third  layer.  The  choking  of  the  flow  at 
two  locations  restricts  the  mass,  momentvim  and  energy  flux  into  the  enclosure.  Thus, 
the  designed  shock-reflector  appears  to  limit  the  high-energy  flow  outside  the  enclosure 
from  entering.  The  pressure  inside  and  outside  the  enclosure  has  to  equalize  eventually. 
However,  the  objective  was  to  design  a  passive  reflector  that  will  slow  this  process  to  a 
rate  that  can  be  tolerated  by  the  sensitive  equipment  inside  the  enclosure.  This  objective 
has  been  achieved. 
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A  final  note  relating  to  the  computational  performance  of  the  code:  computation  time 
for  approximately  12,000  time  steps  consumed  approximately  8  hours  of  CPU  time  on  a 
Cray  YMP  8/128  computer  (single  processor).  Computation  time  per  node  per  time  step 
was  approximately  30  microseconds. 

2.2.1  Shock  Wave  Diffraction  Processes. 

After  examining  the  overall  success  of  this  device  in  mitigating  shock  amplitude  and 
harmonic  content  entering  the  chamber,  we  will  examine,  in  detail,  shock  diffraction  aroimd 
the  first  bar  and  the  first  chevron.  As  noted  above,  the  greatest  advantage  of  the  adaptive 
grid  methodology  is  that  while  it  enables  us  to  obtain  an  engineering  level  solution  (i.e.,  a 
correct  solution  obtained  within  reasonable  CPU  and  memory  expenditures  at  a  fast  turn¬ 
around  pace),  it  also  allows  the  examination  of  the  most  minute  physical  detail,  enabled 
by  the  minimum  mesh  size  of  25  pm. 

2.2.2  Shock  Wave  Diffraction  Around  the  FVont  Bar. 

Figiue  2-9  shows  expanded  views  of  Mach  number,  density  and  vortidty  contours 
around  the  top  of  the  front  bar  between  t=0.1375  ms  to  t=0.21  ms.  To  better  appreciate 
the  resolution  displayed  in  these  figures,  it  is  noted  that  the  bar  thickness  is  1  cm  (10000 
pm).  A  minimum  normal  size  of  25  pm  would  allow  placing  about  400  points  across  the 
bar,  a  sufficient  resolution. 

At  t=0.1375  ms,  about  60  micro-seconds  after  shock  impact,  and  only  about  28  micro¬ 
seconds  after  shock  diffraction  around  the  downstream  comer,  two  fully  roUed-up  vortices 
are  shown  to  be  shed  from  both  comers.  The  vortex  shed  from  the  upstream  comer 
develops  a  classic  Kelvin- Helmholtz  instability,  with  growing  amplitude  in  the  downstream 
direction.  The  multiple  shocklets  shown  at  this  time  (Figure  2-9a  at  t=0.1375)  to  emanate 
from  the  downstream-side  of  this  vortex  resulted  from  the  interaction  of  the  supersonic  flow 
above  the  vortex  with  the  unstable  shear-layer  -  a  flow  resembling  an  vmsteady  sup>ersonic 
flow  over  a  curved  wavy  wall.  It  is  postulated  that  the  system  of  soimd  waves  originated 
from  the  roll-up  process  of  the  shear  layer  (it  is  well  known  that  a  roll-up  vortex  acts  as  a 
quadmpole)  is  amplified  by  the  supersonic  flow.  This  assumption  is  supported  by  the  fact 
that  similar  simulations  with  lower  (subsonic)  Mach  numbers  flows  (and  lower  shock  over¬ 
pressures)  did  not  produce  the  system  of  oriented  (normal  to  the  flow)  shocklets  observed 
here.  This  process  will  be  further  investigated  in  the  future. 

The  two  shed  vortices  continue  to  grow,  almost  self-similarly,  as  shown  by  the  density 
and  vortidty  results  at  t=0.15  ms  (Figures  2-9b  and  2-9g,  respectively).  Eventually,  down¬ 
stream  convection  of  the  sound  produced  by  the  multiple-shocklet  system  resvtlted  in  an 
instability  growth  at  the  downstream  vortex  at  t=0.18  ms,  with  an  identical  frequency,  as 
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shown  by  both  the  density  and  Mach  number  results.  Figures  2-9c  and  2-9d,  respectively. 

The  growing  upstream,  unstable  vortex  reached  the  downstream  comer  at  t=0.21 
ms.  Figures  2-9e,  2-9f  and  2-9h  show  the  Mach  number  and  vorticity  contours,  at  t=0.21 
ms  and  0.245  ms,  respectively.  These  results  show  the  growing  instability  of  the  shear 
layer  originating  from  the  upstream  comer,  the  merging  and  lifting  of  the  two  vortices 
by  the  entrained  flow  below  the  downstream  vortex,  the  strong  shock  attached  to  the 
lifting  upstream  vortex,  the  two  recircxilation  bubbles  within  this  vortex,  and  the  system 
of  shocklets  (or  flow  instability)  still  attached  to  the  bottom  vortex. 

2.2.3  Shock  Wave  Diflfraction  Around  the  First  Chevron. 

The  incident  shock  diffracted  arotmd  the  first  chevron  upstream  comer  at  about 
t=0.25  ms.  The  density  contours  at  t=0.30  ms,  only  50  /is  after  shock  passage  (Fig¬ 
ure  2- 10a),  show  a  fully  developed,  rolled-up  vortex  sheet  attached  to  this  comer.  This 
figure  also  shows  part  of  the  leading  incident  shock,  and  its  reflection  from  the  bottom 
plane.  The  reflected  shock  has  been  broken  by  the  chevron.  Notice  that  the  reflected 
shock  imder  the  chevron  has  propagated  further  than  the  shock  on  the  left,  as  it  has  been 
accelerated  by  the  vortex.  In  addition,  the  reflected  shock  from  both  faces  of  the  chevron 
should  have  been  continuous.  Instead,  on  the  underside  of  the  chevron,  the  reflected  shock 
has  been  broken  by  the  diflerent-direction  velocities  in  the  different  segments  of  the  vortex. 

Figtue  2-lOb  at  t=0.35  ms,  and  Figures  2-lOc  and  2-lOd,  at  t=0.40  ms,  respectively, 
show  the  growth  of  the  vortex  and  the  downstream  propagation  of  the  Kelvin-Helmholtz 
instability  originated  near  the  comer.  The  system  of  imstable  shocks  that  resulted  from 
the  interaction  of  the  mean  supersonic  flow  with  the  convoluted  vortex  attached  to  the 
upstream  bar  (Figure  2-5c  at  t=0.45  ms),  has  reached  this  comer,  as  shown  in  Figures  2- 
11a.  In  addition,  the  primary  shock  reflected  from  the  center  of  the  first  chevron  (Figures 
2-4  and  2-5)  has  diffracted  aroimd  the  comer.  The  results  show  the  break-up  of  this  shock 
into  two  primary  elements:  the  down-propagating  part  outside  the  vortex,  and  the  part 
that  has  diffracted  around  the  comer  into  the  vortex.  This  part  has  been  decelerated  by 
the  vortex  since  it  was  propagating  against  the  flow  near  the  wall.  The  shock  between 
these  cuts  has  been  diffused,  or  sheared,  by  the  local  strong  gradient  (i.e.,  the  subsonic 
to  supersonic  transition)  across  the  vortex  sheet  (Figirre  2-lla2).  In  addition,  another 
shock,  which  has  previously  reflected  from  the  bottom  plane,  is  shown  at  this  time  to  pass 
through  the  vortex,  and  is  located  exactly  at  the  comer.  When  traced  through  the  vortex, 
it  is  shown  that  the  shock  has  been  stretched  or  compressed  depending  on  the  propagation 
direction  relative  to  the  local  flow;  the  shock  immediately  outside  the  vortex  has  been 
accelerated  downstream;  inside  the  vortex  the  shock  was  either  accelerated  or  decelerated, 
depending  on  the  local  mean  flow  direction. 
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Figures  2-1  la  through  2-llg  show  a  series  of  density,  Mach  number  and  vorticity 
contours  at  times  of  0.45,  0.48,  0.49,  0.50,  0.51,  0.55  and  0.60  ms,  respectively.  This 
series  clearly  shows:  a)  continuous  growth  of  the  unstable  vortex;  b)  multiple  shock-vortex 
interactions  that  resulted  in  partial  break-up  of  the  vortex;  and  c)  shock  acceleration- 
deceleration  by  the  local  flow  within  and  outside  the  vortex.  Notice,  for  instance,  the 
pair  of  shocks  that  appeared  on  the  right  side  at  t=0.49  ms.  These  shocks  reflected  from 
the  second  chevron  (Figure  2-5b).  The  ciuvature  of  these  shock  changed  in  time  since 
different  parts  of  the  shock  were  convected  at  different  locaJ  velocities;  the  shock  portion 
near  the  bottom  plane  was  convected  against  mean  flow  Mach  number  significantly  lower 
than  at  the  drctunference  of  the  vortex  (Figure  2-6d).  Another  example  is  shown  in  Figure 
2- 11c  through  g  as  the  shock  that  has  reflected  from  the  bottom  plane  and  propagated 
upward  near  the  edge  of  the  chevron  was  partially  transmitted  and  partially  reflected  by 
the  vortex.  The  portion  outside  had  propagated  upward  slower  than  the  shock  p>ortion 
inside  the  vortex  (Figure  2-llfi),  simply  because  the  mean  flow  outside  was  transonic 
and  in  the  opposite  direction,  while  the  flow  inside  the  vortex  was  almost  stagnant.  Still 
another  example  is  the  second  reflected  shock  from  the  second  chevron,  moving  in  the 
observed  frames  from  right  to  left.  Comparison  of  the  results  at  t=0.55  and  t=0.60  ms 
(Figures  2-llf  and  g,  respectively)  shows  the  stretching  and  compression  of  this  shock  as 
it  passes  through  the  vortex.  While  at  t=:0.55  ms  this  shock  was  almost  continuous  (at 
the  right  side  of  the  figure),  at  t=0.60  ms  the  broken  shock  is  clearly  traceable  through 
the  vortex,  up  to  the  reflection  from  the  bottom  face  of  the  first  chevron;  d)  shock-induced 
vortex  instability.  Shock  interaction  with  the  unstable  vortex  is  shown  to  induce  further 
instability,  and  some  directionality  in  the  vortex  break-up  process.  Notice,  for  instance, 
vortex  break-up  due  to  the  interaction  with  the  upstream  moving  shock  shown  in  Figures 
2-llfi  through  2-llf3,  and  2-llgi  through  2-llg3,  with  the  strongest  effect  immediately 
after  shock  passage  (compare,  for  instance,  Figures  2-llf2  and  2-llg2). 

2.3  SUMMARY  AND  CONCLUSIONS. 

A  recently  developed  transient,  two-dimensional,  finite-element,  shock-capturing  scheme 
on  imstructured  grids  was  applied  to  the  study  of  a  shock  interaction  with  a  passive  shock- 
reflector  placed  at  an  opening  of  a  vented  enclosure.  The  objective  of  this  effort  was  to 
design  a  passive  device  which,  while  allowing  the  normal  ventilation  of  the  enclosure  under 
steady  conditions,  will  prevent  blast-waves  impinging  on  the  wall  from  entering  the  enclo- 
s\ire.  Shock  wave  attenuation  resulted  here  from  the  shock  diffraction  processes  aroimd 
the  chevrons,  that  damped  the  incident  shock  wave  amplitude  and  harmonic  content,  and 
eventually  transformed  the  shock  to  a  low-  amplitude  compression  wave  by  the  time  it 
entered  the  enclosure.  In  addition,  the  sharp  comers  of  the  grill  were  relied  on  to  produce 
vortices  which  significantly  reduced  the  effective  flow  area  (a  self-  choking  device)  to  slow 
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the  outside  quasi-steady  high-pressure  flow  behind  the  shock  from  entering  the  chamber. 

FVom  the  numerical  point  of  view,  the  results  demonstrate  the  capability  of  the  devel¬ 
oped  adaptive  refinement/ coarsening  algorithm  to  properly  adapt  to  weak  shocks,  rarefac¬ 
tion  waves  and  contact  discontinuities  and  the  resultant  excellent  resolution  of  the  captured 
flow  features.  All  shocks  had  propagated  very  long  distances,  but  were  still  captured  as 
sharp,  nonoscillatory  discontinuities.  Grid  refinement  was  observed  as  the  shocks  propa¬ 
gated  into  new  areas,  while  grid  coarsening  was  observed  in  areas  already  traversed  by  the 
shocks.  In  addition,  the  captured  shocks’  thickness  did  not  change  over  several  thousand 
steps,  as  all  shocks  were  still  captured  over  two  to  three  elements,  while  the  element  size 
around  the  shock  has  been  unchanged,  thanks  to  the  adaptive  procedure. 

Among  the  many  interesting  physical  processes  monitored  in  this  computation  were 
the  growth  of  imstable  shear  layers,  transient  shock  interaction  with  unstable  shear  layers, 
shock-shock  interaction,  shock-rarefaction  interaction,  etc.  Vortices  are  observed  to  shed 
from  all  sharp  comers.  The  vortex  roll-up  process  was  enhanced  by  vortex  interaction 
with  the  system  of  curved  reflected  shocks  reverberating  in  the  system.  Since  vortex  dy¬ 
namics  is  a  phenomenon  normally  associated  with  viscous  flow  processes  while  the  present 
numerical  scheme  only  solves  the  Euler  equations,  we  have  investigated  this  phenomenon 
and  concluded  that  vorticity  production  immediately  after  shock  passage  resulted  from  the 
“Bjerknes  Effect.”  The  shed  vortices  resulted  in  effective  flow-area  restrictions.  Specifi¬ 
cally,  the  vortices  shed  from  the  first  and  second  layers  of  chevrons  almost  completely 
blocked  the  flow  between  the  first  and  second  chevrons,  choking  the  flow  and  restricting 
the  mass,  momentum  and  energy  flux  into  the  enclosure.  Thus,  the  designed  shock-reflector 
appears  to  limit  the  high-energy  flow  outside  the  enclosure  from  entering.  The  pressure 
inside  and  outside  the  enclosure  has  to  equalize  eventually.  However,  the  objective  was  to 
design  a  passive  reflector  that  will  slow  down  this  process  to  a  rate  that  cam  be  tolerated 
by  the  sensitive  equipment  inside  the  enclosure.  This  objective  hais  been  achieved. 
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SECTION  3 

NUMERICAL  SIMULATION  OF  SHOCK  INTERACTION 
WITH  COMPLEX-GEOMETRY  3-D  ABOVE-GROUND  STRUCTURES 


This  section  describes  the  application  of  a  recently  developed  three-dimensional  adap¬ 
tive  finite  element  shock  capturing  scheme  on  unstructured  tetrahedral  grids  to  the  sim¬ 
ulation  of  shock  diffraction  processes  around  typical  complex  geometry  three-dimensional 
structures  such  as  tanks  and  missiles.  The  3-D  surfaces  were  defined  using  a  new  CAD- 
CAM-like  user-friendly  solid-body  generator  (PREGEN3D).  The  advancing  front  method 
then  generated  the  volume  grid.  The  shock  diffraction  simulations  were  initiated  by  im¬ 
posing  the  initial  and  boundary  conditions;  only  two  levels  of  refinement  were  used  to 
refine  the  shock  at  its  initial  position.  The  resolution  and  fidelity  of  the  simulated  shock 
wave  diffraction  phenomena,  performed  via  a  solution  of  the  transient  compressible  Eu¬ 
ler  equations,  were  enhanced  by  the  application  of  the  classic  h-enrichment/coarsening 
grid  adaptation  scheme,  with  density  as  the  critical  adaptation  parameter.  A  high  de¬ 
gree  of  vectorization  was  achieved  by  pre-sorting  the  elements  and  then  performing  the 
refinement/coarsening  on  the  assembled  groups.  Further  reductions  in  CPU-requirements 
were  realized  by  optimizing  the  identification  and  sorting  of  elements  for  refinement  and 
deletion. 

The  computational  results  obtained  demonstrate  the  successful  application  of  the  new 
3-D  adaptation  procedure  to  shock  interaction  with  curved  surfaces  -  a  new  capability.  Ex¬ 
cellent  shock  adaptation  and  resolution  are  obtmned  at  all  times;  all  shocks  are  captured 
as  sharp  discontinuities,  without  producing  pre-  or  post-shock  oscillations.  Several  inter¬ 
esting  three-dimensional  shock  diffraction  processes  are  identified  and  discussed  in  detail. 
It  is  shown  that  the  geometry  must  be  described  precisely  in  order  to  obtain  the  correct  lift 
and  drag  forces  on  the  vehicle.  Finally,  the  results  demonstrate  the  robust  p>erformance  of 
the  method  and  show,  at  least  for  the  simulation  of  strongly  imsteady  fiows,  considerable 
savings  in  both  CPU-time  and  storage  over  fixed-mesh  structured  grid  schemes. 

3.1  INTRODUCTION. 

The  solution  of  large-scale  transient  problems  around  complex  geometries  is  a  common 
problem  to  many  areas  of  computational  fiuid  dynamics.  Before  developing  a  numerical 
methodology  to  simulate  these  flows,  one  must  choose  between  structured  or  imstructured 
grid  schemes.  Both  have  their  advantages  and  disadvantages.  The  biggest  disadvantages 
of  the  structured  grid  approach  are:  a)  the  gridding  of  a  reasonably  complex,  engineer¬ 
ing  type  3-D  structure  may  require  many  man-months,  compared  to  days/weeks  with  the 
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unstructured  grid  approach;  and  b)  the  complexity  involved  in  developing  adaptation  algo¬ 
rithms  on  structxired  grids  for  strongly  transient  features.  The  advantages  of  the  structiued 
grid  methodology  over  xmstructured  grids  are:  a)  significantly  reduced  memory  overhead 
requirements  (about  5-15  times  lower);  and  b)  reduced  computation  time  per  node  per 
time  step  (about  3  times  faster).  Thus,  for  a  typical  3-D  computation,  the  unstructiued, 
adaptive  approach  may  be  superior  only  for  flow  simulation  around  a  complex  geometry 
structure  and/or  shock-wave  propagation  over  large  temporal  and  spatial  domains. 

Past  attempts  to  simulate  transient  shock  interaction  with  complex  geometry  3-D  bod¬ 
ies  have  been  limited.  Some  experience  was  gained  applying  a  structured,  non-adaptive 
grid  approach  to  simulate  shock  interaction  with  3-D  semi-submeiged  structures.  One 
approach  prescribed  the  sloped  surfaces  of  the  body  as  ‘‘staircases”  [1,2],  which  would  nat¬ 
urally  yield  a  somewhat  erroneous  solution.  Another  disadvantage  of  this  approach  is  that 
numerical  boundary  layers  are  developed  whenever  shocks  climbing  along  the  staircases 
represent  ctirved  boimdaries. 

The  objective  of  the  present  research  was  to  numerically  simulate  shock  dii&action 
processes  about  complex-geometry  3-D  structures.  Under  a  separate  project,  we  recently 
developed  a  numerical  methodology  capable  of  efficiently  simulating  transient  shock-shock 
and  shock-structme  interactions  for  realistic,  complex,  engineering-type  three-dimensional 
geometries.  The  imstructured  grid  approach  [3]  was  the  most  suitable  for  generating 
surface  and  volume  grids  for  complex-geometry  bodies.  In  addition,  since  our  typical 
areas  of  interest  require  the  simulation  of  strong  shocks,  both  steady  and  transient,  a 
high-resolution  monotonicity-preserving  algorithm  for  unstructured  grids  was  developed. 
The  method,  called  FEM-FCT  [4],  is  based  on  Zalesak's  [5]  generalization  of  the  Flux- 
Corrected  Transport  (FCT)  algorithms  [6,7]  to  multidimensional  problems.  Extensive 
application  of  FEM-FCT  to  2-D  simulations  [8,9],  as  well  as  limited  experience  with  3- 
D  (shock  wave  propagation  over  simplified  geometries  [10]),  has  demonstrated  excellent 
agreement  with  experiment2Ll  data.  With  these  schemes,  both  traveling  as  well  as  stationary 
shocks  are  captured  within  two-three  gridpoints  (for  either  2-D  or  3-D)  without  the  over- 
and  imdershoots  that  appear  in  linear  schemes.  Finally,  since  the  flowfields  typically  are 
smooth  except  for  a  few  regions  where  strong  gradients  appear,  efficient  adaptive  refinement 
techniques  for  transient  problems  are  required. 

As  noted  above,  a  major  advantage  of  the  imstructtired  over  the  structiired  grid  ap¬ 
proach  is  the  ease  with  which  each  complex  geometry  structure  can  be  discretized  [11]. 
Generating  the  surface  or  volume  grids  for  a  typical  airplane  or  vehicle  may  take  several 
man-months  using  a  structured  grid  approach,  and  only  one  to  two  weeks  using  an  imstruc¬ 
tured  grid  [11,12].  A  second  very  attractive  feature  of  unstructvired  grids  is  the  ease  with 
which  adaptive  refinement  can  be  incorporated.  Since  additional  degrees  of  freedom  do 


31 


not  destroy  any  previous  structure,  the  flow  solver  requires  no  further  modiflcation  when 
operating  on  an  adapted  grid.  For  many  practical  problems,  the  regions  that  need  to  be 
reflned  are  extremely  small  as  compared  to  the  overall  domain.  Therefore,  the  savings  in 
storage  and  CPU-requirements  typically  range  between  50-100  as  compared  to  an  overall 
fine  mesh  [8,13,14].  Our  exp>erience  in  2-D  [8,9]  indicates  that  for  the  majority  of  the  daily 
production-type  nms,  adaptive  refinement  makes  the  difference  in  the  ability  or  inability 
to  rtm  the  problems  to  an  acceptable  accuracy  in  a  reasonable  time.  Without  it,  we  would 
be  forced  to  use  much  coarser  grids,  with  lower  accuracy,  for  the  same  expense.  Here  we 
extend  the  2-D  capabilities  to  3-D.  Although  more  complex  in  coding,  the  rewards  of  a 
3-D  adaptive  refinement  capability  are  greater  than  those  encoimtered  in  2-D. 

The  meshes  used  for  the  present  calculations  were  generated  using  FRGEN3D  [12,17]. 
This  imstructured  grid  generator  is  based  on  the  advancing  front  method.  After  defin¬ 
ing  the  siuface  description  of  the  domain  to  be  gridded,  these  siirfaces  are  triangulated. 
Thereafter,  the  face  forming  the  smallest  new  element  is  deleted  from  the  front,  and  a  new 
element  is  added.  This  process  is  repeated  recursively  until  no  more  faces  are  left  in  the 
front.  FRGEN3D  generates  approximately  25,000  tetrahedra  per  minute  on  the  Cray-2. 
This  number  may  increase  for  small  meshes  (less  than  50,000  tetrahedra),  as  more  front 
collapses  per  element  generated.  For  very  large  meshes,  a  global  h-refinement  option  is 
available.  With  it,  the  rate  of  generation  is  increased  to  about  200,000  tetrahedra  per 
minute. 

Over  the  past  three  years  it  became  clear  that  FRGEN3D  by  itself  was  not  sufficient 
to  quickly  generate  the  mesh  required  for  an  arbitrary  3-D  problem.  Although  it  solves  the 
generation  problem  once  and  for  all,  the  surface  generation  then  becomes  the  dominant 
man-hour  bottleneck.  Therefore,  under  a  separate  research  project  we  recently  developed 
FECAD,  a  suite  of  tools  that  allows  the  user  to  produce  FRGEN3D-compatible,  error- 
free  input  in  a  faster  way.  FECAD  not  only  allows  the  user  to  exercise  basic  CAD- 
CAM  operations  (shrinking,  translations,  rotations,  surface  lofting,  etc.),  but  also  eases 
the  merging  of  several  parts  of  the  surfeice  into  one  cohesive,  well-defined  input-file.  This 
allows  the  merger  of  files  produced  by  different  users  and/or  different  surface  generators. 
FECAD  has  a  whole  series  of  built-in  diagnostics  to  avoid  such  undesirable  featwres  as 
doubly  defined  points,  isolated  points  or  lines,  and  badly  defined  lines  or  siirfsu;es.  FECAD 
proved  invaluable  when  trying  to  construct  in  a  matter  of  days  am  error-free  FRGEN3D- 
compatible  input-file. 

Another  important  capability  that  was  developed  over  the  past  four  years  was  PRE- 
BACK,  a  semi-structured  3-D  background  generator.  Experience  with  the  user  community 
has  shown  that  the  preferred  way  of  generating  background  grids  is  by  lofting  a  2-D  tri- 
angidation  in  the  third  dimension  or  by  rotation.  PREBACK  allows  both  operations.  In 
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addition,  it  allows  interactive  operations  to  move  background  grid-points  aroimd  (trans¬ 
lation,  rotation,  shrinking,  etc.),  as  well  as  to  modify  the  grid-generation  parameters  in 
space  (size,  shape). 

The  flow  solver  employed  was  FEFL074,  a  3-D  ALE  hydro-solver  based  on  FEM-FCT 
[10,15,16].  H-reflnement  [13,16]  is  the  preferred  approach  for  grid  adaption.  The  high  order 
scheme  used  is  the  consistent-mass  Taylor-Galerkin  algorithm.  Combined  with  a  modifled 
second  order  Lapidus  artiflcial  viscosity  scheme,  the  resiilting  scheme  is  second  order  accu¬ 
rate  in  space  and  fourth  order  accurate  in  phase.  The  spatio-temporal  adaptation  is  based 
on  local  h-refinement,  where  the  reflnement/deletion  criterion  is  a  modified  H2-seminoim 
[13].  Based  on  past  experience  with  simulations  of  shock  wave  propagation  processes  in 
both  2-D  and  3-D,  density  was  chosen  as  the  critical  parameter  for  the  refinement /deletion 
criterion. 

Post-processing  was  performed  with  the  FEPOST3D  and  MOVIESUBS  packages. 
FEPOST3D  is  based  on  FEPLOT4D  [18],  but  performs  all  the  CPU-intensive  filtering 
operations  on  the  Cray.  Only  the  plane-  or  surface-triangulations  me  sent  back  to  the 
SGI-IRIS-4D/80GT  for  plotting.  The  user  specifies  before  the  run  the  planes  and  surfaces 
to  be  inspected.  Although  seemingly  trivial,  this  step  was  in  a  large  part  responsible  for 
the  success  of  the  present  effort.  It  implied  waiting  30  seconds  for  a  plot-file  (about  2.5- 
4Mbytes  of  information),  instead  of  hours  (about  130-160Mbytes  for  a  complete  flowfield 
mesh).  In  addition,  it  allowed  the  production  of  movies. 

In  the  following  sections  we  will  briefly  review  some  of  the  main  schemes  that  are  in¬ 
tegrated  into  FEFL074,  such  as  the  efiicient  error  indicator  and  the  FCT-FEM  algorithm. 
More  detailed  analysis  is  foimd  in  Ref  [19].  In  addition,  we  include  a  description  of  some 
of  the  recent  algorithm  modifications  and  improvements  that  reduced  CPU  requirements 
and  improved  the  performance  of  the  methodology. 

3.2  THE  ERROR  INDICATOR. 

Many  variants  of  an  efficient  error  indicator  had  been  investigated  in  the  past  [3- 
10].  Here,  we  seek  a  method  that  is  efficient  and  reliable  for  transient  compresrible  flow 
problems.  This  leads  us  to  the  following  design  criteria  for  the  error  indicator: 

a)  The  error  indicator  should  be  fast. 

b)  The  error  indicator  should  be  dimensionless,  so  that  several  key  variables  can  be 
monitored  simultaneously. 

c)  The  error  indicator  should  be  boimded,  so  that  no  further  user  intervention  becomes 
necessary  as  the  solution  evolves. 
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d)  The  error  indicator  should  not  only  mark  the  re^ons  with  strong  shocks  to  be  refined, 
but  also  weak  shocks,  contact  discontinuities  and  other  weak  features  in  the  flow. 

For  the  refinement  method,  the  design  criteria  are  as  follows: 

a)  The  method  should  be  conservative,  i.e.,  a  mesh  change  shovild  not  result  in  the 
production  or  loss  of  mass,  momentum  or  energy. 

b)  The  method  should  not  produce  elements  that  are  too  small,  as  this  woiild  reduce  too 
severely  the  allowable  timestep  of  the  explicit  flow  solvers  employed. 

c)  The  method  shotild  be  fast.  In  particvdar,  it  should  lend  itself  to  some  degree  of 
parallelism. 

d)  The  method  should  not  involve  a  major  storage  overhead. 

3.2.1  The  Error  Indicator. 

An  error  indicator  that  meets  the  design  criteria  (a-d)  was  proposed  in  [3]  as  follows: 


error  s= 


\steoni  itTivaiivts\ 
h  {first  derivatives]  +  c  |meon  vo/tte| 


(3.1) 


By  dividing  the  second  derivatives  by  the  absolute  value  of  the  first  derivatives,  the  error 
indicator  becomes  boimded,  dimensionless,  and  the  “eating  up”  effect  of  strong  shocks  is 
avoided.  The  terms  following  c  are  added  as  a  noise  filter  in  order  not  to  refine  “wiggles” 
or  “ripples”  which  may  appear  due  to  loss  of  monotonicity.  The  value  for  e  thus  depends 
on  the  algorithm  chosen  to  solve  the  PDEs  describing  the  physical  process  at  hand.  The 
multidimensional  form  of  this  error  indicator  is  given  by 


= 


(3.2) 


where  denotes  the  shape-function  of  node  I.  This  error  indicator  has  performed  very 
well  in  2-D  over  the  years  [3,4].  However,  when  first  used  in  3-D,  it  proved  unreliable.  The 
sotirce  for  this  seemingly  inconsistent  behavior  was  found  to  originate  from  the  large  local 
variations  in  element  size,  shape,  as  well  as  the  niunber  of  elements  surrounding  a  point. 
These  will  produce  large  variations  of  the  second  term  in  the  denominator  that  are  not 
based  on  physics,  but  on  the  mesh  structiure  itself.  The  solution  was  to  modify  the  error 
indicator  as  follows: 
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(3.3) 


E'  = 


E.  ,(/o  |W' +  €MMAT,kT^IU,\' 


where  MMATj  is  the  lumped  mass-matrix  at  point  I,  and  hj  the  average  element  length 
at  point  I.  This  error  indicator  proved  to  be  remarkably  insensitive  to  local  variations 
in  element  size  and  shape,  while  still  yielding  the  correct  indicator  values  for  physical 
phenomena  of  interest.  We  attribute  this  robust  performance  to  the  smoothing  effects 
of  two  averaging  operations  working  simultaneotisly:  the  liunped  mass-matrix  and  the 
point-lengths. 


After  determining  the  values  of  the  error  indicators  in  the  elements,  all  elements 
l}ring  above  a  preset  threshold  value  CTORE  are  refined,  while  all  elements  lying  below  a 
pr  threshold  value  CTODE  are  coarsened.  Protective  layers  of  elements  are  added  to 
SI  id  the  elements  to  be  refined,  so  that  the  “feature”  (e.g.,  a  shock)  always  travels  into 

an  already  refined  region.  The  number  of  protective  layers  that  are  added  depends  on  the 
Courant-number  employed  and  the  number  of  time  steps  taken  between  grid  modifications. 
Usually  the  refinement  is  performed  every  5-10  time  steps,  so  that  a  Courant-number  of 
Cn=0.2  is  sufficient. 


3.2.2  Adaptive  Refinement  Method. 

Our  previous  experience  in  2-D  indicates  that  the  only  two  refinement  methods  that 
are  tnily  general  and  efficient  for  the  class  of  problems  considered  here  are  h-refinement 
[3,5]  and  remeshing  [7,8,10].  However,  for  strongly  unsteady  problems,  where  a  new  grid  is 
required  every  5-10  time  steps,  local  h-refinement  seems  to  be  preferable.  Several  reasons 
can  be  given  for  this  choice: 

a)  Conservation  presents  no  problem  for  h-refinement. 

b)  No  interpolations  other  than  the  ones  nattirally  given  by  the  element  shape-functions 
are  required.  Therefore,  no  numerical  diffusion  is  introduced  by  the  adaptive  refine¬ 
ment  procedure.  This  is  in  contrast  to  adaptive  remeshing,  where  the  grids  before  and 
after  a  mesh  change  may  not  have  the  srune  points  in  common.  The  required  inter¬ 
polations  of  the  unknowns  will  result  in  an  increased  amount  of  numerical  diffusion 

19]. 

c)  H-refinement  is  very  well  suited  to  vector-  and  parallel  processors.  This  is  of  particular 
importance  in  the  present  context,  where  a  mesh  change  is  performed  every  seven  time 
steps,  and  a  large  percentage  of  mesh  points  is  affected  in  each  mesh  change. 

d)  H-refinement  is  more  robust  than  remeshing. 
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As  described  above,  we  limit  the  number  of  refinement/coarsening  levels  per  mesh  change 
to  one.  Moreover,  we  only  allow  refinement  of  a  tetrahedron  into  two  (along  a  side), 
four  (along  a  face)  or  eight  new  tetrahedra.  We  call  these  tetrahedra  1:2,  1:4  and  1:8 
tetrahedra  or  refinement  cases  respectively.  At  the  same  time,  a  1:2  or  1:4  tetrahedron  cam 
only  be  refined  further  to  a  1:4  tetrahedron,  or  by  first  going  back  to  a  1:8  tetrahedron 
with  subsequent  further  refinement  of  the  8  sub-elements.  We  call  these  the  2:4,  2:8-t-  and 
4:8+  refinement  cases.  The  refinement  cases  are  summarized  in  Figure  3-1.  This  restrictive 
set  of  refinement  rules  is  necessary  to  avoid  the  appearance  of  ill-deformed  elements.  At 
the  same  time,  it  considerably  simplifies  the  refinement /coarsening  logic.  An  interesting 
phenomenon  that  does  not  appear  in  2-D  is  the  apparently  free  choice  of  the  inner  diagonal 
for  the  1:8  refinement  case.  As  shown  in  Figure  3-2,  we  can  place  the  inner  four  elements 
aroimd  the  inner  diagonals  5-10, 6-8,  or  7-9.  In  the  present  case,  the  shortest  inner  diagonal 
was  chosen.  This  choice  produces  the  smallest  amoimt  of  distorted  tetrahedra  in  the  refined 
grid.  When  coarsening,  we  again  only  allow  a  limited  number  of  cases  that  are  compatible 
with  the  refinement.  Thus,  the  coarsening  cases  become  8:4,  8:2,  8:1,  4:2,  4:1,  2:1.  These 
coarsening  cases  are  summarized  in  Figure  3-3. 

When  constructing  the  algorithm  to  refine  or  coarsen  the  grid,  one  faces  the  usual 
decision  of  speed  versus  storage.  The  more  information  from  the  previous  grid  that  is 
stored,  the  faster  the  new  grid  may  be  constructed.  As  storage  requirement  minimization 
was  one  of  the  goals  of  the  present  research,  we  tried  to  keep  only  the  essential  information 
needed  between  mesh  changes  while  minimizing  CPU  time.  This  was  accomplished  by  using 
a  modified  tree-structure  which  requires  twelve  integer  locations  jjer  element  in  order  to 
identify  the  parent  and  son  elements  of  any  element,  as  well  as  the  element  type. 

The  first  seven  integers  store  the  new  elements  (sons)  of  an  element  that  has  been 
subdivided  into  eight  (1:8).  For  the  1:4  and  1:2  cases,  the  sons  are  also  stored  in  this 
allocated  space,  and  the  remaining  integer  locations  are  set  to  zero.  In  the  eighth  integer 
the  element  from  which  the  present  element  originated  (the  parent  element)  is  stored.  The 
ninth  integer  denotes  the  position  number  in  the  parent  element  from  which  this  element 
came.  The  tenth  integer  denotes  the  element  type.  We  can  either  have  parents  or  sons  of 
1:8, 1:4  or  1:2  tetrahedra.  We  mark  these  with  a  positive  value  of  the  element  type  for  the 
parents,  and  a  negative  value  for  the  sons.  Thus,  for  example,  the  son  of  a  1:8  element 
would  be  marked  as  -8.  Fintdly,  in  the  eleventh  and  twelfth  integer  locations  the  local  and 
global  refinement  levels  are  remembered.  These  twelve  integer  locations  per  element  are 
sufficient  to  construct  further  refinements  or  to  reconstruct  the  original  grid. 

3.2.3  Algorithmic  Implementation. 

Having  outlined  the  basic  refinement /coarsening  strategy,  we  can  now  describe  in  more 
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depth  its  algorithinic  implementation.  One  complete  grid  change  algorithmically  requires 
the  following  five  steps: 

1)  Construction  of  the  missing  grid  information  needed  for  a  mesh  change. 

2)  Identification  of  the  elements  to  be  refined. 

3)  Identification  of  the  elements  to  be  deleted. 

4)  Refinement  of  the  grid  where  needed. 

5)  Coarsening  of  the  grid  where  needed. 

1.  Construction  of  Missing  Grid  Tnfornmtinn.  The  missing  information  consists  of  the 
sides  of  the  mesh  and  the  sides  belonging  to  each  element.  The  sides  are  dynamically 
stored  in  two  arrays,  one  containing  the  two  points  each  side  connects  and  the  other  one  (a 
pointer-array)  containing  the  lowest  side-number  reaching  out  of  a  point.  The  formation  of 
these  two  arrays  is  accomplished  in  three  miun  loops  over  the  elements,  which  are  partially 
vectorizabk.  After  having  formed  these  two  side-arrays,  a  further  loop  over  the  elements 
is  performed,  identifying  which  sides  belong  to  each  element. 

2.  Identification  of  Elements  to  be  Refined.  The  aim  of  this  sub-step  is  to  determine  on 
which  sides  further  gridpoints  need  to  be  introduced,  so  that  the  resulting  refinement 
patterns  on  an  element-level  belong  to  the  allowed  cases  listed  above,  thus  producing  a 
compatible,  valid  new  mesh.  Five  main  steps  are  necessary  to  achieve  this  goal: 

a)  Mark  elements  that  require  refinement; 

b)  Add  protective  layers  of  elements  to  be  r^ned; 

c)  Avoid  elements  that  become  too  small,  or  that  have  been  refined  too  often. 

d)  Obtain  preliminary  list  of  sides  where  new  points  will  be  introduced; 

e)  Add  further  sides  to  this  list  imtil  an  admissible  refinement  pattern  is  achieved. 

a)  Mark  elements  that  require  refinement 

Using  the  modified  error  indicator  given  by  E)q.  (1.3)  and  the  prescribed  refinement  toler¬ 
ance  CTORE,  those  elements  that  require  further  refinement  are  marked  using  the  array 
LELEM(1:NELEM).  This  array  is  marked  as  LELEM(IE)=1  ^  element  is  to  be  refined, 
LELEM(IE)=0  element  is  not  to  be  refined. 

b)  Add  protective  layers  of  elements  to  be  refined 

If  protective  layers  of  elements  are  to  be  added  ahead  of  the  featture  to  be  refined,  we 
perform,  for  each  additional  layer,  the  following  set  of  operations: 
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-  Zero  an  integer  point-array  (e.g.,  LPOIN(l;NPOIN)=Oy, 

-  Loop  over  the  elements  to  be  refined,  marking  (e.g.,  LP0IN(IP)=1)  all  points  con¬ 
nected  to  these  elements; 

-  Zero  the  integer  element-array  (e.g.,  LELEM(1;NELEM)=0)-, 

-  Loop  over  all  elements;  if  at  least  one  point  of  a  given  element  has  been  marked,  refine 
this  element  (e.g.,  LELEM(IE)=1) 

c)  Avoid  elements  that  become  too  small,  or  that  have  been  refined  too  often 

A  sharp  feature  in  the  flow  domain,  like  a  shock,  will  tend  to  produce  error  indicator 
values  that  always  lie  above  the  refinement  tolerance  CTORE.  As  a  consequeace,  elements 
close  to  sudi  a  feature  will  be  refined  every  time  the  mesh  is  adapted.  In  order  to  avoid 
this  “refinement  ad  infinitum,"  one  has  to  impose  either  a  maximum  permissible  number 
of  refinement  levels  per  element,  or /and  a  minimum  allowable  element  size.  Given  these 
constraints,  those  elements  which  are  already  too  small  (if  a  minimum  allowed  element 
size  has  been  given),  or  have  already  been  refined  too  many  times  (if  a  maximum  allowed 
nmnber  of  refinement  levels  has  been  prescribed),  are  deleted  firom  the  list  of  elements  to 
be  refined. 

d)  Obtain  preliminary  list  of  sides  for  new  points 

Given  the  side/element  information  obtained  above,  we  can  now  determine  a  first  set  of 
sides  on  which  new  gridpoints  need  to  be  introduced.  This  set  of  sides  is  still  preliminary, 
as  we  only  allow  certain  ty|>es  of  refinement. 

e)  Add  further  sides  to  this  list  until  an  admissible  refinement  nattem  is  achieved 

The  list  of  sides  marked  for  the  introduction  of  new  points  is  still  preliminary  at  tins 
point.  In  most  cases,  it  will  not  lead  to  an  admissible  refinement  pattern  to  construct  a 
new  mesh.  Therefore,  further  sides  are  marked  for  the  introduction  of  new  points  until  an 
admissible  refinement  pattern  is  reached.  This  is  accomplished  by  looping  several  times 
over  the  elements,  checking  on  an  element  level  whether  the  set  of  sides  marked  can  lead 
to  an  admissible  new  set  of  sub-elements.  The  algorithm  used  is  based  on  the  observation 
that  the  admissible  cases  are  based  on  the  introduction  of  new  points  along  one  side  (1:2), 
three  contiguous  sides  (1:4),  or  six  contiguous  sides  (1:8).  These  admisable  cases  can  be 
obtained  from  the  following  element-by-element  algorithm: 

-  Set  the  node-array  LNODE(1:4)=0; 

-  Loop  over  the  sides  of  the  element:  if  the  side  has  been  marked  for  the  introduction  of 
a  new  point,  set  LN0DE(IP1)=1,  LN0DE(IP2)=1,  where  IPl,  IP2  are  the  end-nodes 
corresponding  to  this  side; 

-  Loop  over  the  sides  of  the  element:  if  LN0DE(IP1)=1  and  LN0DE(IP2)=1,  mark 
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the  side  marked  for  the  introduction  of  a  new  point. 

Practical  calculations  with  several  admissible  layers  of  refinement  and  large  grids  revealed 
that  sometimes  up  to  15  passes  over  the  mesh  were  required  to  obtain  an  admissible  set  of 
sides.  A  considerable  reduction  in  CPU  was  achieved  by  presorting  the  elements  as  follows: 

-  Add  up  all  the  sides  marked  for  refinement  in  an  element; 

-  If  0,1  or  6  sides  were  marked:  do  not  consider  further; 

-  If  4  or  5  sides  were  marked:  mark  all  sides  of  this  element  to  be  refined; 

-  If  2  or  3  sides  were  marked:  analyze  in  depth. 

This  then  yields  the  final  set  of  sides  on  which  new  gridpoints  are  introduced. 

3.  Identification  of  Elements  to  be  Deleted.  The  mm  of  this  sub-step  is  to  determine  which 
points  are  to  be  deleted,  so  that  the  resulting  coarsening  patterns  on  an  element-level  belong 
to  the  allowed  cases  listed  above,  thus  producing  a  compatible,  valid  new  mesh.  Four  main 
steps  are  necessary  to  achieve  this  goal: 

a)  Mark  elements  to  be  deleted; 

b)  Filter  out  elements  where  the  parent  and  all  sons  are  to  be  deleted; 

c)  Obtain  preliminary  list  of  points  to  be  deleted; 

d)  Delete  points  from  this  list  imtil  an  admissible  coarsening  pattern  is  achieved. 

a)  Mark  elements  to  be  deleted 

As  before,  we  start  by  determining  from  the  modified  error  indicator  given  by  Elq.  (1.3)  and 
the  prescribed  deletion  tolerance  CTODE,  those  elements  that  should  be  coarsened.  Thus, 
we  mark  an  element  array  LELEM(1:NELEM)  as  follows:  LELEM(IE)=-1  element  is 
to  be  deleted,  LELEM(IE)=0  =»  element  is  not  to  be  deleted. 

b)  Filter  out  elements  where  parent  and  all  sons  are  to  be  deleted 

It  is  clear  that  only  those  elements  should  be  deleted,  for  which  both  the  parent  as  well 
as  all  its  sons  have  been  marked  for  deletion.  Therefore,  only  the  parent  elements  to  be 
coarsened  are  considered  further.  For  these  elements,  a  check  is  performed  whether  their 
resp>ective  son  elements  have  also  been  marked  for  deletion.  If  any  of  the  son  elements  have 
not  been  marked  for  deletion,  neither  the  parent-element  nor  any  of  its  sons  are  considered 
further. 

c)  Obtain  preliminary  list  of  points  to  be  deleted 

Given  the  list  of  parent-elements  to  be  coarsened,  we  can  now  determine  a  preliminary 
list  of  points  to  be  deleted.  Thus,  all  the  points  that  would  be  deleted  if  all  the  elements 
contained  in  this  list  were  coarsened  are  marked  as  ‘‘total  deletion  p>oints”. 
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d)  Delete  points  from  this  list  until  an  admissible  coarsening  pattern  is  achieved 

The  list  of  total  deletion  points  obtained  in  the  previous  step  is  only  preliminary,  as  imal- 
lowed  coarsening  cases  may  appear  on  an  element  level.  We  therefore  perform  loops  over 
the  elements,  deleting  all  those  total  deletion  points  that  would  result  in  allowed  coarsen¬ 
ing  cases  for  the  elements  adjoining  them.  This  process  is  stopjjed  when  no  incompatible 
total  deletion  points  are  left. 

4.  Refinement  of  the  Grid  Where  Needed.  The  introduction  of  further  points  and  elements 
is  performed  in  two  independent  steps,  which  in  principle  could  be  performed  in  parallel. 

1)  Points:  To  add  further  points,  the  sides  marked  for  refinement  in  sub-step  2  are  grouped 
together.  For  each  of  these  sides  a  new  grid-point  will  be  introduced.  The  interpolation  of 
the  coordinates  and  unknowns  is  then  performed  using  the  side/point  information  obtained 
in  sub-step  1.  These  new  coordinates  and  imknowns  are  added  to  their  respective  arrays. 
In  the  same  way  new  boimdary  conditions  are  introduced  where  required. 

2)  Elements:  In  order  to  add  further  elements,  the  sides  marked  for  refinement  are  labelled 
with  their  new  gridpoint-number.  Thereafter,  the  element/side  information  obtained  in 
sub-step  1  above  is  employed  to  add  the  new  elements.  The  elements  to  be  refined  are 
grouped  together  according  to  the  refinement  cases  shown  in  Figme  3-1.  Each  case  is 
treated  in  block  fashion  in  a  separate  subroutine.  Perhaps  the  major  breakthrough  of  the 
present  work  was  the  reduction  of  the  many  possible  refinement  cases  to  only  six.  In  order 
to  accomplish  this,  some  information  for  the  2:8+  and  the  4:8+  cases  is  stored  ahead  in 
scratch  arrays.  After  these  elements  have  been  refined  according  to  the  2:8  and  4:8  cases, 
their  sons  are  screened  for  further  refinement  using  this  information.  All  sons  that  require 
further  refinement  are  then  group>ed  together  as  1:2  or  1:4  cases,  and  processed  in  turn. 

5.  Coarsening  of  the  Grid  Where  Needed.  The  deletion  of  points  and  elements  is  again 
p>erformed  in  two  indei>endent  steps,  which  in  principle  could  be  performed  in  parallel. 

a)  Points:  The  points  to  be  deleted  having  been  marked  in  sub-step  3  above,  all  that 
remains  to  be  done  is  to  fill  up  the  voids  in  the  coordinate-,  unknown-  and  boimdary 
condition-arrays  by  renumbering  points  and  boundary  conditions. 

b)  Elements:  The  deletion  of  elements  is  again  performed  blockwise,  by  grouping  together 
all  elements  corresponding  to  the  coarsening  cases  shown  in  Figure  3-3.  Thereafter,  the 
elements  are  also  renumbered  (in  order  to  fill  up  the  gaps  left  by  the  deleted  elements), 
and  the  point-renumbering  is  taken  into  consideration  within  the  connectivity-arrays. 

3.2.4  Recent  Algorithm  Modifications. 

We  have  recently  significantly  modified  the  adaptive  finite  element  methodology  for 
transient  problems  in  3-D  [6].  Since  for  typical  shock  calculations,  mesh  adaption  takes 
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place  every  5-10  time  steps,  we  concluded  that  every  stage  of  the  adaption  process  must  be 
thoroughly  optimized.  The  two  areas  that  required  the  most  intense  optimization  efforts 
were: 

a)  Error  indicators  for  grids  with  large  local  variations  of  element  size  and  shape,  and 

b)  Faster  construction  of  the  new  mesh. 

Improved  error  indicators:  As  an  error  indicator  we  use  a  modified  interpolation  theory 
error  indicator: 
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where  /  lenotes  the  shape-function  of  node  I.  By  dividing  the  second  derivatives  by  the 
absolute  value  of  the  first  derivatives  the  error  indicator  becomes  bounded,  dimensionless, 
and  the  ‘eating  up’  effect  of  strong  shocks  is  avoided.  This  error  indicator  has  performed 
very  well  in  2-D  over  the  years  [2].  However,  when  first  used  in  3-D,  it  proved  iinreliable. 
The  source  for  this  seemingly  inconsistent  behavior  was  found  to  stem  from  the  large  local 
variations  in  element  size,  shape,  as  well  as  the  number  of  elements  surroimding  a  point 
encouintered  in  typical  3-D  unstructured  grids.  These  will  produce  large  variations  of  the 
second  term  in  the  denominator,  which  are  not  based  on  physics  but  on  the  mesh  structure 
itself.  The  solution  was  to  modify  this  error  indicator  as  follows: 
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where  MMATj  is  the  lumped  mass-matrix  at  point  /,  and  hj  the  average  element  length 
at  point  /.  This  error  indicator  proved  to  be  remarkably  insensitive  to  local  variations 
in  element  size  and  shape,  while  still  yielding  the  correct  indicator  values  for  physical 
phenomena  of  interest.  We  attribute  this  good  performance  to  the  smoothing  effects  of 
two  averaging  operations  working  simultaneously:  the  lumped  mass-matrix  imd  the  point- 
lengths. 

Faster  construction  of  the  new  mesh. 

The  main  CPU-intensive  operations  performed  during  one  mesh  change  are  finding  the 
sides  and  the  faces  of  each  element,  determining  the  refinement  and  coarsening  patterns, 
correcting  boundary  p>oints,  and  renumbering  the  elements.  Although  seemingly  trivial, 
these  op)erations  account  for  a  significant  p>ercentage  of  total  CPU  time.  We  developed 
new,  optimal  algorithms  for  them.  Two  examples  are  described  here  in  more  depth: 
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1)  Determining  the  list  of  sides  for  new  points 

Given  the  side/element  information,  and  a  list  of  elements  to  be  refined,  a  first  set  of 
sides  on  which  new  gridpoints  need  to  be  introduced  is  determined.  In  most  cases,  it 
will  not  lead  to  an  admissible  refinement  pattern.  Therefore,  further  sides  are  marked 
for  the  introduction  of  new  points  until  an  admissible  refinement  pattern  is  reached. 
This  is  done  by  looping  several  times  over  the  elements,  checking  on  an  element  level 
whether  the  set  of  sides  marked  can  lead  to  an  admissible  new  set  of  sub-elements. 
Practical  calculations  revealed  that  sometimes  up  to  15  passes  over  the  mesh  where 
required  to  obtain  an  admissible  set  of  sides.  An  80%-90%  reduction  in  CPU  was 
achieved  by  presorting  the  elements  as  follows; 

-  Add  up  all  the  sides  marked  for  refinement  in  an  element; 

-  If  0,1  or  6  sides  were  marked:  do  not  consider  further; 

-  If  4  or  5  sides  were  marked:  mark  all  sides  of  this  element  to  be  refined; 

-  If  2  or  3  sides  were  marked:  analyze  in  depth. 

2)  Element  renumbering 

In  order  to  vectorize  the  element  assembly  as  much  as  possible,  the  elements  are 
renumbered,  such  that  within  each  assembly  pass  a  point  is  accessed  only  once  by 
the  elements.  This  renumbering  has  to  take  place  after  every  mesh  change.  Before 
optimization,  it  took  over  10%  of  the  total  CPU  time  for  typical  runs.  The  renumbering 
subroutine  was  optimized  by  working  only  on  the  remaining  elements,  and  extensive 
scalar  optimization,  minimizing  the  number  of  operations  and  memory  access.  With 
the  new  renumbering  algorithm,  the  CPU-time  required  for  this  operation  dropp>ed  to 
less  than  1%. 

3.2.5  Recent  Modifications  to  the  H-Refinement  Algorithm. 

Over  the  past  two  years,  we  have  significantly  modified  the  adaptive  finite  element 
methodology  for  transient  problems  in  3-D  [10].  Since  for  typical  shock  calculations  mesh 
adaption  takes  place  every  5  to  10  time  steps,  it  has  been  concluded  that  every  stage  of 
the  adaption  process  must  be  thoroughly  optimized.  The  two  areas  that  required  the  most 
intense  optimization  efforts  were  error  indicators  for  grids  with  large  local  variations  of 
element  size  and  shape,  and  faster  construction  of  the  new  mesh. 

Improved  error  indicators. 

As  an  error  indicator  we  use  a  modified  interpolation  theory  error  indicator; 
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where  denotes  the  shape-function  of  node  I.  By  dividing  the  second  derivatives  by  the 
absolute  value  of  the  first  derivatives,  the  error  indicator  becomes  bounded,  dimensionless, 
and  the  ‘^eating  up”  effect  of  strong  shocks  is  avoided.  This  error  indicator  has  performed 
very  well  in  2-D  over  the  years  [8].  However,  when  first  used  in  3-D,  it  proved  unreliable. 
The  source  of  this  seemingly  inconsistent  behavior  was  found  to  stem  from  the  large  local 
variations  in  element  size  and  shape,  as  well  as  the  number  of  elements  surrounding  a  point 
encountered  in  typical  3-D  unstructured  grids.  These  wil  produce  large  variations  of  the 
second  term  in  the  denominator  that  are  not  based  on  physics,  but  on  the  mesh  structure 
itself.  The  solution  was  to  modify  this  error  indicator  as  follows: 
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(3.7) 


where  MMATj  is  the  lumped  mass-matrix  at  point  I,  and  hj  the  average  element  length 
at  point  I.  This  error  indicator  proved  to  be  remarkably  insensitive  to  local  variations 
in  element  size  and  shape,  while  still  yielding  the  correct  indicator  values  for  physical 
phenomena  of  interest.  We  attribute  this  good  performance  to  the  smoothing  effects  of 
two  averaging  operations  working  simultaneously:  the  lumped  mass-matrix  and  the  point- 
lengths. 

Faster  construction  of  the  new  mesh. 

The  main  CPU-intensive  operations  performed  during  one  mesh  change  are  finding  the 
sides  and  the  faces  of  each  element,  determining  the  refinement  and  coarsening  patterns, 
correcting  boundary  jrnints,  and  remunbering  the  elements.  Although  seemingly  trivial, 
these  operations  accoimt  for  a  significant  percentage  of  total  CPU  time.  We  developed 
new,  optimal  algorithms  for  them.  Two  examples  are  described  here  in  more  depth: 

a)  Determining  the  List  of  Sides  for  New  Points 


Given  the  side/elem^t  information,  and  a  list  of  elements  to  be  refined,  a  first  set  of 
sides  on  which  new  gridpoints  need  to  be  introduced  is  determined.  In  most  cases,  it 
will  not  lead  to  an  admissible  refinement  pattern.  Therefore,  further  sides  are  marked 
for  the  introduction  of  new  points  until  an  admissible  refinement  pattern  is  reached. 
This  is  done  by  looping  several  times  over  the  elements,  checking  on  an  element  level 
whether  the  set  of  sides  marked  can  lead  to  an  admissible  new  set  of  sub-elements. 
Practical  calcrilations  revealed  that  sometimes  up  to  15  passes  over  the  mesh  were 
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required  to  obtain  an  admisable  set  of  sides.  An  80%-90%  reduction  in  CPU  was 
achieved  by  presorting  the  elements  as  follows: 

-  Add  up  all  the  sides  marked  for  refinement  in  an  element; 

-  If  0,1  or  6  sides  were  marked;  do  not  consider  further; 

-  If  4  or  5  sides  were  marked;  mark  all  sides  of  this  element  to  be  refined; 

-  If  2  or  3  sides  were  marked;  analyze  in  depth, 
b)  Element  Renumbering 

In  order  to  vectorize  the  element  assembly  as  much  as  possible,  the  elements  are 
renumbered  such  that  within  each  assembly  pass  a  point  is  accessed  only  once  by 
the  elements.  This  renvimbering  has  to  take  place  after  every  mesh  change.  Before 
optimization,  it  took  over  10%  of  the  total  CPU  time  for  typical  runs.  The  renumbering 
subroutine  was  optimized  by  working  only  on  the  remaining  elements,  and  extensive 
scalar  optimization  by  minimizing  the  number  of  operations  and  memory  access.  With 
the  new  renumbering  algorithm,  the  CPU-time  required  for  this  operation  drof^d  to 
less  than  1%. 

3.3  SHOCK  INTERACTION  WITH  A  GENERIC  TANK. 

The  first  simulation  of  shock  diffraction  about  a  complex-geometry  three-dimensional 
structxire  conducted  under  this  project  was  a  simulation  of  shock  interaction  with  a  generic 
tank  configtiration  (closely  resembling  the  West  German  Leopard).  Since  the  shock  im¬ 
pinged  head-on  on  the  tank,  a  plane  of  symmetry  exists  and  only  half  of  the  domain  can 
be  modelled. 

Initial  tests  with  the  new  three-dimensional  code  described  above  demonstrated  its 
ability  to  capture  moving  and  stationary  shocks  over  three  elements  without  loss  of  mono¬ 
tonicity.  For  this  class  of  problems  and  the  algorithm  employed  it  was  foimd  that  the 
following  choice  of  refinement/coarsening  parameters  produced  nonoscillatory  shocks  and 
rarefactions: 

-  refinement  tolerance:  CT0RE=Q.\S 

-  coarsening  tolerance:  CTODE^SSXQ 

-  noise  parameter:  c  =0.12 

-  key  variable:  density 

-  refinement  frequency:  every  7  time  steps 

-  number  of  protective  layers:  none 

-  Courant  number  of  the  hydro-solver:  C  =  0.45. 
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Figures  3-4  through  3-12  show  the  adapted  computational  grid  and  pressure  contours 
at  several  times  during  the  computation.  Since  it  is  difficult  to  visualize  results  in  the 
computational  volume,  results  will  be  presented  only  on  the  boundaries  or  on  specified 
planes.  The  -l-X  direction,  the  shock  propagation  direction  (Figure  3-4),  will  be  referred 
to  as  the  axial  direction,  the  -l-Y  direction  as  the  vertical  direction,  and  the  -4-Z  direction 
as  the  cross-stream  direction.  Plane  X-Y  is  the  plane  of  symmetry,  as  shown  in  Figure  3-4. 
Since  it  was  intended  to  integrate  the  solution  to  t=50  ms,  the  computational  domain  had 
to  be  significantly  larger  than  shown  in  this  figiure:  20  m  upstream  of  the  tank  (to  allow 
for  shock  reflection),  25  m  downstream  of  the  tank,  15  m  in  the  vertical  direction  and  20  m 
in  the  cross-stream  direction.  Past  experience  with  2-D  simulations  [4]  demonstrated  that 
advantage  of  the  adaptive,  unstructxued  grid  approach  vs.  the  structured  grid  approach 
when  applied  to  shock  propagation  over  large  distances.  The  ability  to  adapt  and  refine 
only  where  needed  was  shown  to  be  even  more  important  for  the  3-D  simulations.  The 
shocks  developed  in  the  computational  domain  are  always  captured  as  sharp  discontinuities, 
demonstrating  the  effectiveness  of  the  new  methodology.  The  results  (es|)ecially  at  later 
times)  demonstrate  the  advantage  of  iising  an  adaptive  scheme,  as  the  reflected  3-D  shocks 
were  sharply  captiured  even  after  propagating  significant  distances.  A  similarly  resolved 
shock  computed  with  a  fixed,  structtired  mesh  would  have  required  maintaining  a  fine  grid 
resolution  everywhere  in  the  computational  volume,  at  a  prohibitively  high  computational 
cost.  Grid  refinement  is  observed  as  the  shocks  propagate  into  new  volumes,  while  grid 
coarsening  is  observed  in  areas  already  traversed  by  the  shocks.  As  expected,  the  later 
time  results  exhibit  significantly  weaker  shocks  than  observed  in  the  past  for  corresponding 
2-D  simulations,  due  to  three-dimensional  expansion  processes. 

Figure  3-4  shows  a  superimposition  of  the  pressure  contoxus  and  the  adapted  mesh 
on  the  boundaries  at  t=0.  A  shock  with  an  amplitude  of  10  psi  overpressure  and  low 
supersonic  Mach  number  (1.14)  was  placed  approximately  1.0  m  upstream  of  the  cannon. 
The  initial  mesh  consisted  of  71,219  points  and  405,400  elements.  Later  in  the  computation, 
the  average  computational  mesh  consisted  of  approximately  230,000  points  and  1.3  million 
tetrahedra. 

The  shock  impinged  on  the  cannon  at  t=:3.0  ms  (Figure  3-5),  at  which  time  a  basically 
planar  reflection  was  observed.  Three-dimensional  expansion  processes  would  immediately 
relieve  the  high  pressure  reflected  shock.  Thus,  at  t=4.2  ms  (Figinre  3-6b),  the  three- 
dimensionally  expanded  shock  was  observed  at  a  X-Y  plane  located  0.5  m  from  the  plane 
of  symmetry,  while  the  amplitude  of  the  reflected  shock  has  been  drastically  reduced  and 
was  considerably  lower  than  predicted  by  2-D  models.  The  computational  mesh,  demon¬ 
strated  excellent  adaptation  to  the  3-D  shock  (Figures  3-6a  and  3-7a,  at  t=4.2  and  5.5 
ms,  respectively).  Three  levels  of  refinement  are  observed,  with  a  minimum  normal  height 
of  0.6  cm.  The  shock  is  about  to  impinge  on  the  front  chassis  of  the  tank  at  t=5.5  ms. 
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Nevertheless,  the  large  tetrahedra  ahead  of  the  shock  (Figure  3-7a),  between  the  shock  and 
the  treads,  were  refined  in  preparation  for  shock  arrival,  and  thus  the  refined  mesh  on  the 
front  chassis  and  treads.  The  pressure  contours  at  this  time  (Figure  3-7b)  show  further 
expansion  of  the  refiected  shock  from  the  cannon,  although  the  refiected  shock  has  not  yet 
arrived  to  a  X-Y  plane  located  1.5  m  from  the  plane  of  symmetry. 

The  propagating  shock  hits  the  front  chassis  at  t=5.9  ms.  The  adapted  mesh  and 
pressure  contours  on  the  tank  surface,  shown  in  Figures  3-8a  and  3-8b,  respectively,  demon¬ 
strate  the  grid  adaptation  to  fiow  gradients.  In  addition  to  the  well-refined  incident  and 
refiected  shocks,  strong  expansions  are  observed  around  all  sharp  comers.  The  incident 
shock  is  shown  on  the  chassis  and  the  sides  of  the  tank.  Several  3-D  reflected  shocks  are 
observed.  Figures  3-8c  and  3-8d  show  pressure  contoiurs  on  three  X-Y  planes,  located  at 
Z=0.5  m,  Z=1.5  m  and  Z=2.75  m,  respectively.  The  results  at  all  planes  exhibit  an  almost 
spherical  reflection  (shock  1)  from  the  slanted  front  of  the  chassis  (shown  in  Figure  3-8a). 
The  reflected  shock  has  impacted  on  the  ground.  Figure  3-8e  shows  a  superimposition 
of  the  pressure  contours  on  the  boundaries  and  two  X-Y  planes  located  at  Z=1.5  m  and 
Z=2.75  m  from  the  plane  of  symmetry  (shown  in  the  figure  on  the  left  and  right  sides, 
respectively).  Although  the  incident  shock  was  planar  and  the  reflection  in  the  vertical 
direction  was  fairly  uniform  above  the  chassis,  the  reflection  pattern  below  the  chassis 
is  more  complex.  While  the  incident  shock  passed  under  the  tank  unimpeded  (Figiires 
3-8c  and  3-8d),  strong  reflections  from  the  treads  are  observed  (Figure  3-8e).  The  higher 
stagnation  pressures  ahead  of  the  treads  resulted  in  greater  upstream  reflection  distances 
than  near  the  plane  of  symmetry,  where  the  incident  shock  weis  free  to  propagate  under 
the  tank  and  the  observed  upstream  reflection  is  from  the  chassis.  It  should  be  noted 
that  the  outline  of  the  shocks  shown  in  Figure  3-8e  represent  the  outer  surface  of  a  single, 
very  complex  three-dimensional  reflected  shock  (called  shock  1).  New  graphic  tools  are 
needed  to  exhibit  the  three-dimensional  natvure  of  these  shocks.  We  intend  to  develop  such 
capabilities  in  the  near  future. 

The  reflected  shocks  from  the  treads  me  three-dimensional  and  later  expEuid  to  fill  the 
space  vmder  the  chassis.  Thus,  at  later  times,  a  more  uniform  (i.e.,  planar)  reflection  front 
was  observed  to  propagate  upstream  (Figure  3-lOb).  The  reflection  of  shock  1  from  the 
ground  (shock  2)  impinged  on  the  treads  and  reflected  (shock  3  in  Figure  3-8b).  While  the 
incident  shock  propagated  past  the  sloped-down  panel  of  the  front  chassis,  a  downward 
propagating  shock  reflected  and  tailed  the  incident  shock  (shock  D  in  Figure  3-8c).  A 
strong  expansion  arotmd  this  comer  is  shown  in  Figures  3-8c  and  3-8d.  Examination 
of  these  figures  (note  Figures  3-8b  and  3-8e)  indicates  that  the  reflected  shock  from  the 
chassis  (shock  1),  with  a  fairly  high  pressvue  amplitude,  impinged  on  the  cannon,  exerting 
additional  bending  moments. 
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The  incident  shock  impacted  on  the  turret  at  about  t=:11.9  ms  (Figure  3-9a).  Strong 
reflection  firom  the  turret  (shock  Ti)  is  observed  at  t=12.87  (Figure  3-9b).  Since  the 
rarefaction  waves  propagated  from  the  turret  edges  to  the  center,  the  largest  pressure 
amplitude,  and  hence  furthest  upstream  reflection  distance,  was  observed  on  the  chassis  at 
the  plane  of  symmetry  (Fig\ire  3-9b).  Notice  that  at  this  time  the  incident  shock  continued 
forward  propagation  in  the  narrow  spacing  between  the  turret  and  the  chassis,  and  has 
almost  impacted  on  the  cylindrical  turret  base. 

The  complex  system  of  three-dimensional  shocks  shown  in  Figures  3-8  continued  their 
evolution.  The  primary  reflected  shock  (shock  1)  continued  expanding,  while  the  pressure 
rise  across  the  shock  continued  to  decrease.  The  decay  rate  was  significantly  higher  than 
for  corresponding  2-D  calculations,  no  doubt  due  to  the  added  expansion  around  the  tank 
in  the  cross-stream  (Z)  direction.  While  shock  3  continued  expanding  ahead  of  the  treads 
(Figure  3-9b),  shock  2  (tmder  the  chassis)  reflected  from  the  chassis  (shock  2a  in  Figure 
3-9c),  and  shock  D  reflected  from  the  ground,  propagated  up  and  was  about  to  impact  on 
the  bottom  of  the  tank. 

The  incident  shock  reflected  from  the  front  of  the  cylindrical  turret  base  at  t=13.0  ms. 
The  adapted  mesh  and  presstire  contours  figures  at  t=15.53  ms,  Figures  3-lOa  and  3- 10b, 
respectively,  demonstrate  continued  mesh  adaptation  to  all  flow  gradients,  shocks  as  well  as 
rarefactions.  Figures  3- 10b  and  3-11  (and  many  others  not  shown  here)  demonstrate  strong 
reflection  (shock  T2)  from  the  turret  base  at  the  plane  of  symmetry.  The  shock  amplitude 
and  reflection  distance  decreased  with  increased  angle  from  the  plane  of  symmetry,  as 
predicted  by  previous  simxilations  of  planar  shock  interaction  with  a  cylinder  [4].  Similarly, 
the  reflection  from  the  turret  (shock  Ti)  expanded  spherically.  Figure  3-11  (at  t=15.53  ms) 
attempts  to  present  the  spherical  expansion  of  shock  Ti,  as  demonstrated  by  the  pressure 
contoiurs  on  the  plane  of  symmetry  (an  X-Y  plane),  the  timk  chassis  (an  X-Z  plane),  and  a 
Y-Z  plane  at  2.0  m  (shock  Tu).  Meanwhile,  complex  shock  evolution  continued  imder  the 
tank,  as  shown  in  Figure  3- 10c,  for  a  X-Y  plane  at  Z=0.5  m.  While  the  primary  reflected 
shock  (shock  1)  continued  upstream  reflection,  shock  2  expanded  in  the  vertical  direction. 
Its  reflection,  combined  with  the  reflection  of  shock  D,  from  the  down-slop>ed  siuface  at 
the  front  of  the  tank,  is  shock  2b.  The  reflection  of  the  curved  shock  2a  (Figure  3-9c)  from 
the  ground  is  shock  2c,  while  shock  2d  is  the  reflection  of  shock  D  from  the  bottom  of  the 
tank  and  the  groimd. 

The  calculation  was  continued  to  later  times.  Due  to  the  complexity  of  the  evolved 
solution,  the  latest  results  presented  are  for  tsl9.69  ms,  when  the  shock  was  still  propagat¬ 
ing  over  the  turret  (Figure  3-12a).  The  back-view  (Figure  3-12b)  shows  a  nice  expansion 
of  the  shock  behind  the  circular  turret  base,  and  demonstrates  the  sharp,  nonoscillatory 
incident  and  reflected  shocks. 
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Total  computation  time  (to  t=25  ms  real  time)  was  approximately  23-25  CPU  hours 
on  a  128  MWords  Cray  2  computer.  Computation  time  per  node  per  time  step  was 
approximately  180  microseconds,  including  the  refinement/deletion  process,  which  added 
approximately  20-25  percent  to  the  overall  computation  time. 

Finally,  it  is  noted  that  a  deficiency  of  the  xmstructured  grid  schemes  is  the  large 
memory  overhead  required  for  the  three-dimensional  algorithm:  approximately  380  words 
per  node.  As  noted  above,  this  overhead  can  be  reduced  to  less  than  250  words  per  node 
by  re-computing  certain  geometric  variables  rather  than  storing  them.  Nevertheless,  such 
an  overhead  necessitates  the  usage  of  the  large  static-memory,  relatively  slow,  Cray  2  (128 
MWords)  supercomputer  for  most  practical  3-D  simulations.  To  alleviate  this  problem, 
a  new  dynamic  domain-splitting  type  algorithm  will  be  developed  in  the  near  future.  In 
this  algorithm,  the  computational  domain  will  be  split  at  each  refinement/deletion  cycle 
to  multiple  domains.  This  will  allow  the  usage  of  faster,  smaller  static-memory  computers 
with  a  solid-state-disc,  such  as  the  Cray  Y-MP. 

3.4  SHOCK  INTERACTION  WITH  A  T-62  TANK. 

The  second  simulation  conducted  vmder  this  project  was  the  simulation  of  3-D  shock 
diffraction  phenomena  aroimd  a  modem  main  battlefield  tank,  in  this  case  a  T-62  tank.  For 
this  class  of  problems  and  the  algorithm  employed,  past  experience  (and  experience  gained 
during  this  study)  indicated  that  the  following  chcnce  of  refinement /coarsening  parameters 
produced  the  best  results: 

-  refinement  tolerance:  CTORE=Q.QQb 

-  coarsening  tolerance:  CTODE—Q.QI&lh 

-  noise  parameter:  c  =  0.28 

-  key  variable:  density 

-  refinement  frequency:  every  5-7  time  steps 

-  number  of  protective  layers:  none 

-  Courant  number  of  the  hydro-solver:  C  =  0.8. 

The  results  of  the  two  simulations  conducted  in  which  a  strong  shock  impacted  either 
back-on  or  side-on  a  T-62  modem  main  battlefield  tank,  are  shown  in  Figmes  3-13  -  3-19. 
The  pressure  contour  results  are  shown  on  the  surfaces  only  in  the  immediate  vicinity  of 
the  tank,  in  order  to  reduce  the  dxunp-files  size  (dumped  every  10  time  steps),  and  hence 
reduce  field  transfer  and  post-processing  time. 

Figures  3-13a-d  show  several  views  of  the  solid  body  model  and  the  surface  triangula¬ 
tion.  While  the  turret  data  were  available  from  a  CAD-CAM  system  (BRL-CAD),  all  other 
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data  were  measiired  directly  from  the  tank.  The  generation  of  the  surface  and  volume  grids 
took  three  days,  demonstrating  the  versatility  of  the  grid  generation  package  (FRGEN3D, 
PREGEN3D,  PREBACK).  It  should  be  noted  that  once  the  surface  triangulation  is  per¬ 
formed,  the  generation  of  the  volume  mesh,  using  the  advsincing  front  algorithm  [12,17]  is 
completely  automatic.  Generating  such  a  grid  with  a  structured  grid  approach  would  have 
taken  significantly  longer. 

3.5  BACK  IMPACT. 

The  first  simulation  modeled  shock  impact  on  the  tank  from  the  back.  Taking  advan¬ 
tage  of  the  symmetry  of  this  problem,  the  X-Y  plane  was  defined  as  a  plane  of  symmetry. 
The  initial  grid  included  only  20,591  points,  100,343  elements,  and  7,816  surface  points. 
The  incident  shock  was  placed  (t=0)  1.5  meters  behind  the  tank.  After  adapting  to  the 
initiei  shock  using  only  two  levels  of  refinement,  the  number  of  points  increased  to  41,846, 
the  number  of  elements  to  221,482  and  the  number  of  surface  points  to  9,159.  During 
the  computation,  the  number  of  grid  points  increased  to  about  285,000  (the  maximum  we 
could  fit  in  a  Cray-2  with  an  allowable  memory  usage  of  100  million  words),  the  number 
of  elements  to  about  1.8  million,  and  the  number  of  stuface  points  to  about  40,000. 

Supersonic  boundary  conditions  were  imposed  on  all  outflow  boundaries  due  to  the 
high  Mach  number  of  the  flow  behind  the  shock.  Flow  conditions  on  the  upstream  boimdary 
were  prescribed.  Symmetric  conditions  were  posed  at  the  plane  of  symmetry.  Tangential 
boimdary  conditions  were  imposed  on  the  tank.  The  computational  domain  stretched  a 
distance  of  100  m  in  the  axial  direction  (+X),  25  m  in  the  transverse  direction  (+Z), 
and  25  m  in  the  vertical  direction  (+Y).  The  minimum  tetrahedron  normal  size  allowed 
after  two  refinements  was  0.5  cm.  Maintaining  this  resolution  with  a  structured,  fixed- 
mesh  algorithm  on  such  a  large  computational  volume,  even  for  a  multi-zone  (multi-block) 
approach,  would  have  required  a  significantly  larger  (100  times  larger)  computational  grid, 
thus  significantly  increasing  computational  costs  and  memory  requirements. 

Shock  evolution  aroimd  the  tank  during  the  diffraction  phase  is  shown  in  Figures 
3-14a  through  i.  Figure  3-14a  shows  pressure  contours  on  the  tank’s  surface  and  the 
plane  of  symmetry  after  200  time  steps.  Figures  3- 15a  and  b  show  pressure  contours  at 
this  time  on  two  X-Y  planes  located  at  Z=0.5  and  1.5  meters,  respectively.  The  results 
show  the  almost-planar  reflection  of  the  shock  from  the  tail  of  the  back-deck.  Very  little 
pressure  relief  through  the  expansion  of  rarefaction  waves  emanating  from  the  comers  is 
observed,  due  to  the  very  short  time  after  shock  impact  (Figure  3-15a).  Thus,  the  reflected 
pressure  amplitude  at  the  center  of  the  back-plane  was  in  excellent  agreement  with  the 
corresponding  value  obtained  via  the  solution  of  the  1-D  Riemann  problem.  Strong  planar 
shock  reflection  is  also  observed  from  the  back  of  the  ammunition  train.  A  significantly 
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weaker  spherical  reflection  was  observed  from  the  curved  back  of  the  treads  (Figure  3-15b). 
Shock  propagation  over  the  back-deck  shows  no  increase  in  pressure  (from  the  incident 
pressure  value)  as  the  incident  shock  is  normal  to  the  back-deck  and  passed  the  sharp  90" 
back-deck  comer  without  producing  a  Mach-stem. 

Figure  3- 14b  at  400  steps  shows  a  minor  imaged  view  of  pressure  contours  on  the 
surface  of  the  tank.  Pressure  contours  on  the  ground  show  the  reflected  shock  from  the 
back-plane  (still  very  weak  due  to  the  mostly  planar  reflection),  spherical  shock  reflection 
from  the  back  of  the  treads,  and  shock  wave  focusing  between  the  treads  and  the  groimd. 
This  focusing  effect  produced  the  highest  overpressure  values  obtained  in  this  computation. 
This  figure  also  exhibits  one  of  the  more  interesting  processes  observed  in  this  computa¬ 
tion:  shock  diffraction  and  focusing  aroimd  and  between  the  wheels.  High  overpressure 
values  were  obtained  whenever  the  shock  wave,  which  propagated  normal  to  the  primary 
axis  of  the  wheels,  was  focused  between  the  treads  and  a  wheel.  Hence  a  system  of  rever¬ 
berating  shocks  between  the  wheels  and  between  the  wheels/treads/wheel  cover  assembly 
was  formed.  Similarly,  high  overpressure  values  were  obtained  as  the  shocks  converged 
within  the  downstream  arcs  of  the  cylindrical  reams.  Conversely,  low  overpressures  woe 
obtained  within  the  upstream  arcs  of  the  reams,  contributing  to  a  high  torque  force  on  the 
wheels.  Figures  3- 14c  and  3- 15c  show  this  phenomenon  at  t=600  steps.  These  results  also 
show  the  (by  now)  almost  cylindrical  shock  reflection  from  the  back-plane  with  a  large 
pressure  imprint  on  the  groimd  (Figure  3- 15c).  We  observe  here  the  initiation  of  a  Mach 
stem,  with  the  stem,  the  original  reflected  shock,  its  reflection  from  the  ground  and  the 
resulting  triple  point,  all  clearly  defined. 

The  incident  shock  impacted  on  the  turret  at  t=800  steps  (Figure  3-14d).  The  planar 
data  at  Z=1.0  m  (Figure  3-15d)  demonstrate  the  rise  (off  the  ground)  of  the  triple  point 
and  the  truly  three-dimensional  character  of  the  Mach  stem,  as  the  pressure  is  observed 
to  decay  behind  the  stem  with  increased  distance  in  the  transverse  direction.  Analysis 
of  pressure  data  at  other  planes  showed  that  the  overpressure  value  behind  the  stem  was 
fairly  uniform  in  the  +Z  direction  up  to  Z=1.2  m,  about  the  width  of  the  back-deck,  and 
decayed  rapidly  afterwards.  This  conclusion  was  supported  by  later  results  (Figure  3-14e 
at  1000  steps,  and  Figures  3-15e-g  at  1200  steps).  Another  interesting  3-D  shock  reflection 
phenomenon  was  the  reflection  from  the  treads-giound  comer.  The  high  pressure  relief  is 
again  truly  three-dimensional,  as  seen  in  the  pressure  contours  data  on  the  ground,  the 
treads,  and  the  different  planar  data  (Figtues  3-15c  through  3-15f  at  t=1200  steps). 

Complex  shock  diffraction  processes  were  observed  due  to  the  shock  interaction  with 
the  turret.  Figure  3-14e  at  1000  steps,  and  Figxires  3-14g,  3-15e,  f  and  g  at  1200  steps  and 
later  results  demonstrate  that  large  amplitude  reflections  were  obtained  on  turret  surfaces 
normal  or  almost  normal  to  the  flow,  while  weak  reflections  and  rarefactions  were  observed 
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at  turret  surfaces  parallel  or  inclined  at  small  angles  to  the  flow.  The  initial  reflection 
from  the  turret  looked  like  a  bow  shock  at  and  near  the  plane  of  symmetry  (Figmes  3-14e 
and  3-15e),  while  no  reflection  was  observed  near  the  edge  (in  the  +Z  direction)  of  the 
turret,  where  the  incident  shock  was  parallel  to  the  turret  surfsM:e  (Figures  3-15g  and  h  at 
Z=1.0  m  and  1.5  m,  respectively).  At  later  times  (Figures  3-14f  through  h),  pressure  relief 
near  the  top  inhibited  the  upstream  propagation  of  the  reflected  shock,  while  the  shock 
inunediately  above  the  back-deck  propagated  a  much  larger  distsmce  upstream,  forming 
a  very  complex  3-D  shock  with  a  strong  spatial  dependence.  Indeed,  pressure  data  at 
points  on  top  of  the  back-deck  located  only  about  40-60  cm  apart  experienced  significantly 
different  reflected  shock  amplitudes,  as  will  be  shown  later. 

These  results  demonstrate  that  a  precise  geometric  description  of  the  complex  geome¬ 
try  3-D  body  is  required  to  obtain  the  correct  load  on  the  structure.  A  simplified  simulation 
of  this  process  using  a  block  to  represent  the  turret  would  yield  erroneous  loading  on  the 
deck  top:  a  critical  requirement. 

The  primary  reflected  shock  from  the  back-deck  has  reached  an  apparent  “quasi¬ 
steady”  stand-off  distance  after  about  2,200  time  steps.  This  stand-off  distance  is  critically 
dependent  on  the  exact  geometric  description  of  the  deck-treads-side  traines  assembly  that 
determines  the  percentage  of  blocked  area.  The  stand-off  distance,  in  tiun,  determines  the 
lift  and  drag  forces  on  the  tank  during  the  diffraction  and  drag  phases. 

The  results  at  later  times  (after  the  shock  propagated  past  the  turret)  are  shown  in 
Figure  3-14i  (2,800  time  steps).  These  results  (and  others  not  shown  here)  indicate  that 
the  shock  that  emerged  from  under  the  deck  (up-front),  led  the  shock  that  propagated 
above  and  around  the  turret  by  as  much  as  one  meter.  This  higher  propagation  velocity 
and  stagnation  pressure  imder  the  tank  resulted  from  the  propagation  of  the  incident 
shock  across  a  sloped  siirface  at  the  bottom  back  of  the  deck,  a  sinface  that  acted  as  a 
(2-D)  compression  ramp.  The  results  also  show  the  complex  quasi-steady  reflected  3-D 
shock  upstream  of  the  turret.  Flow  expansion  downstream  of  the  tank’s  back-plane  has 
accelerated  the  flow  to  transonic  speeds  near  the  back-deck,  and  to  supersonic  speeds 
high  above  the  deck.  Thus  the  observed  compression  wave  near  the  bottom  of  the  turret, 
and  the  oblique  shock  near  the  top.  Naturally,  this  phenomena  varied  in  the  transverse 
direction  due  to  the  reduced  blockage  by  the  turret. 

Figures  3- 16a  through  d  show  superpositions  of  the  pressure  contours  and  adpated 
mesh  after  200,  400,  600,  and  1000  time  steps,  respectively.  These  results  demonstrate 
the  excellent  shock  adaption  capability  of  the  new  algorithm,  as  observed  on  all  surfaces. 
Investigation  of  several  planar  surfaces  demonstrated  equally  successfiil  adaption  in  the 
voliune.  Only  two  levels  of  refinement  were  utilized  here.  All  shocks  are  observed  to  be 
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fully  refined  and  are  captured  over  two  to  three  elements.  No  wiggles  are  observed  in  the 
solution. 

The  results  presented  here  show  the  shock  evolution  arotmd  the  tank  only  during  the 
diffraction  phase.  The  computation  wras  actually  continued  to  a  later  time,  when  the  shock 
has  propagated  a  long  distance  (about  30  meters)  downstream.  Even  after  propagating 
such  a  long  distance,  all  shocks,  reflected  or  transmitted,  were  still  captured  as  sharp, 
nonoscillatory  discontinuities,  with  a  minimum  tetrahedron  normal  size  of  less  than  1.5 
cm.  These  resvilts  demonstrate  the  advantage  of  the  adaptive,  imstructured  grid  approach 
over  fixed-grid,  structured  codes  (i.e.,  significant  reductions  of  CPU  time  and  memory 
requirements).  Utilizing  these  schemes  to  propagate  a  shock  over  such  long  distances  with 
the  presecribed  resolution  may  be  prohibitively  expensive. 

The  computation  of  more  than  7,000  time  steps  consumed  approximately  60  CPU 
hours  on  a  Cray-2. 

3.6  SIDE  IMPACT. 

The  third  computation  simulated  shock  impingement  side-on  on  the  same  tank.  Fig¬ 
ures  3- 17a  through  c  show  the  superimposed  mesh  and  pressure  contours  at  0,  300,  and  600 
time  steps.  Excellent  grid  adaption  was  obtained.  The  Gourard  shaded  pressure  contours 
at  300,  600,  and  800  time  steps  (Figures  3-18a,  c  and  d,  respectively),  and  the  enlarged 
view  of  the  pressiue  contours  results  between  wheels  three  and  four  at  300  steps  (Figure 
3- 18b)  demonstrate,  in  addition  to  the  high  fidelity  of  the  computations,  many  interesting 
shock  difiraction  phenomena.  Among  them: 

a)  The  stand-off  distance  of  the  quasi-steady  reflected  shock  was  significantly  smaller  than 
for  the  back-on  impingement.  This  reduction  results  &om  the  l2u:ge  spacing  between  the 
wheels  that  allows  the  transmission  of  a  significant  portion  of  the  incident  shock  energy. 
This  phenomenon  points  out  the  importance  of  the  accurate  description  of  the  complex 
geometry  wheels- treads-deck-cover  assembly.  Improp>er  geometric  modeling  that,  for  in¬ 
stance,  will  over-restrict  the  flow  (i.e.,  smaller  spacing),  will  result  in  higher  percentage  of 
reflected  energy,  and  hence  higher  reflected  pressTire  amplitude;  this,  in  turn,  will  resvdt  in 
higher  drag  force  exerted  on  the  vehicle. 

b)  Higher  stagnation  pressure  under  the  deck  than  above  the  turret  resulted  in  higher 
propagation  speed  below,  and  hence,  generated  a  lift  force  on  the  tank. 

c)  The  integrated  drag  force  on  the  tturet  was  higher  thzm  for  the  back-on  attack.  The 
combined  higher  lift  and  drag  forces  will  increase  the  probability  of  the  tank  overttiming. 
Nevertheless,  the  load  exerted  on  the  top  of  the  back-deck  was  significantly  diminished,  as 
the  reflected  shock  amplitude  decreased;  very  little  reflection  from  the  turret  was  observed 
on  the  back-deck  near  the  centerplane  after  800  steps. 
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d)  Due  to  shock  reflection  from  the  turret  and  its  impingement  on  the  cannon,  a  Isurge 
bending  moment  was  exerted  on  the  cannon. 

e)  FHill  stagnation  was  obtained  in  front  of  the  blocked  wheels  (Figures  3-18c  and  d).  The 
shock  traversed  the  outside  of  the  wheels  at  grazing  angles  and  was  only  partially  stagnated 
on  the  side  of  the  main  deck. 

Pressmre-time  histories  at  several  stations  were  compared  to  experimental  data  ob¬ 
tained  in  a  field  test.  The  agreement  between  the  measured  and  predicted  data  was  very 
good. 

The  computation  of  the  diffraction  phase  for  the  side-on  impingement  consumed  ap- 
prcodmately  16  CPU  hours  on  a  Cray-2. 

3.7  COMPARISON  OF  PRESSURE-TIME  HISTORIES. 

Comparisons  of  pressure-time  histories  between  the  results  obtained  for  the  back-and 
side-impact  simulations  are  shown  in  Figiires  3-19a-i  for  several  locations  around  the  tank. 
All  pressure  and  time  data  shown  are  nondimensional.  Figure  3- 19a  shows  pressure-time 
histories  at  the  back-end  of  the  tank.  The  initial  back-impact  shock  completely  stagnated 
to  a  value  that  was  in  very  good  agreement  with  both  1-D  analytical  analysis  and  the 
experimental  data.  The  small  pressure  increase  at  about  t=5  is  the  reflected  wave  from 
the  turret.  In  comparison,  the  shock  for  the  side  impact  traversed  the  back-plate  at  a 
grazing  angle  and  no  reflection  or  pressure  increase  was  observed. 

The  next  four  figures  show  pressure-time  histories  at  foxir  stations  on  the  back-deck. 
The  pressures-contours  data  (Figures  3-14  and  3-15)  demonstrate  the  large  spatial  pressure 
variation  on  the  back-deck  due  to  the  reflected  shock  from  the  turret.  Figures  3-19b  and  c 
show  pressure-time  histories  at  two  stations  on  the  centerline,  one  (station  29)  fairly  close 
to  the  turret,  the  other  (station  34)  only  0.8  m  further  upstream.  For  the  back-impact 
shock,  the  reflected  shock  amplitude  was  reduced  by  about  a  factor  of  three  over  this 
distance  (18.2  at  station  29  vs.  6.5  at  station  34).  For  the  side-impact,  the  differences  are 
not  as  significant,  and  in  fact,  the  station  further  from  the  turret  demonstrates  a  second 
shock  that  resulted  from  shock  diffraction  through  the  spacing  between  the  ammunition 
traines  and  the  back-deck.  Stations  off  the  centerline  demonstrate  shock  overpressure 
dependence  on  distance  from  both  the  turret  and  the  centerline.  Pressure  data  at  station 
43,  located  30  cm  upstream  of  station  34,  and  50  cm  off  centerline,  demonstrate  a  significant 
reduction  in  reflected  shock  overpresstu'e;  the  peak  value  was  about  4.5  times  lower  than 
that  obtained  at  a  station  located  on  the  centerline  at  identical  distance  from  the  turret. 
Pressure  resiilts  at  station  57  (Figure  3-19e),  located  about  1.0  m  further  upstream  and  20 
cm  closer  to  C.L.,  show  further  reduction  of  the  reflected  shock  amplitude.  The  results  for 
the  side-impact  show  a  second  compression  resulting  from  shock  diffraction  through  the 
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spacing  between  the  ammunition  traines  and  the  back-deck.  These  results,  combined  with 
the  pressure-contour  results  (Figures  3-14  and  3-15),  demonstrate  the  need  to  accurately 
model  the  turret.  Modeling  of  the  turret  as  a  square  block  would  have  resulted  in  erroneous 
load  distribution  on  the  back-deck,  and  hence,  erroneous  damage  assessment. 

Figure  3-19f  shows  pressure- time  histories  at  a  station  located  close  to  the  comman¬ 
der’s  hatch  (on  the  turret)  on  a  surface  almost  normal  to  the  side-impact  shock,  but  about 
parallel  to  the  back-impact  shock.  The  geometry  hence  dictates  a  reflected  shock  with  a 
fairly  large  amplitude  for  the  side-impact  case,  in  contrast  to  an  almost  grazing  angle  and 
zero  reflection  for  the  back-impact  shock. 

Pressure  data  at  station  10  (Figure  3-19g),  located  on  the  front  forward-sloping  portion 
of  the  turret  on  the  centerline,  show  expansion  of  the  shock  coming  from  the  back,  with  a 
lower  pressure  than  the  incident  pressure.  The  data  for  the  side-impact  show  that  although 
the  station  was  located  on  a  plane  almost  parallel  to  the  flow,  shock  diffraction  around  the 
rounded  turret  increased  the  pressure  due  to  the  formation  of  a  Mach  stem.  Figure  3-19h, 
for  station  3,  located  on  the  down-sloped  front  deck  just  below  the  cannon,  shows  a  fairly 
significant  shock  reflection  from  the  cannon  for  the  side-impact  case.  The  results  for  the 
back-impact  only  show  the  arrival  of  the  incident,  expanding  shock. 

Finally,  pressure-time  histories  for  station  84  (Figure  3-19i),  located  under  the  tank 
on  the  centerline  at  about  two-thirds  the  length  from  the  front,  show  (for  the  back-impact 
case)  the  arrival  of  the  increased-amplitude  shock,  and  a  smooth  decay  afterwards.  The 
shock  amplitude  increased  due  to  the  compression  by  the  sloped  ix)rtion  of  the  back- 
deck  and  the  formation  of  a  Mach  stem  as  the  shock  propagated  along  the  bottom  of 
the  tank.  For  the  side-imp>act  the  solution  is  even  more  complex.  Since  this  station  was 
located  behi’id  a  wheel,  it  was  not  directly  influenced  by  the  incident  shock  but  rather  by 
the  diffracted  shocks  between  the  wheels.  The  predicted  oscillations  that  followed  were 
the  diffraction  of  the  incident  shock  entering  the  cavity  under  the  deck  through  openings 
between  other  wheels.  Finally,  the  large  pressure  increase  at  about  t=8  was  the  shock 
reflection  from  the  downstream  wheels. 

3.8  SHOCK  INTERACTION  WITH  A  MISSILE  IN  FLIGHT. 

The  last  3-D  shock  difliraction  computation  conducted  under  this  project  simulated 
shock  interaction  at  an  angle  of  45°  with  a  misale  in  flight  (at  an  angle  of  2°).  The 
results  are  shown  in  Figure  3-18.  Here  we  note  that:  a)  due  to  the  3-D  nature  of  the 
shock  diflraction  process,  strong  3-D  rarefaction  effects  are  observed  on  the  surface  of 
the  missile;  b)  the  slope  discontinuities  between  the  stages  produced  rarefaction  waves; 
c)  strong  overpressure  relief  was  observed  at  angles  higher  than  30°;  d)  high  overpressure 
values  were  limited  to  the  cone  zone;  and  e)  excellent  shock  adaption  and  resolution  were 
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obtained  at  all  times;  all  shocks  were  captured  as  sharp  discontinuities,  without  producing 
pre-  or  post-shock  oscillations. 

3.9  SUMMARY  AND  CONCLUSIONS. 

This  section  describes  experience  gained  in  applying  a  recently  developed  three- 
dimensional  adaptive  finite  element  shock  capturing  scheme  on  unstructured  tetrahedral 
grids,  to  the  simulation  of  shock  wave  diffraction  phenomena  aroimd  complex-geometry 
three-dimensional  structures.  The  3-D  surfaces  were  defined  using  a  new  CAD-CAM-like 
user-friendly  solid  body  generator  (PREGEN3D)  that  includes  a  new  automatic  background- 
grid  generator.  The  volumetric  grid  generation  (tetrahedra)  from  the  surface  triangulation 
was  accomplished  using  the  advancing  front  method,  and  was  rated  at  about  25,000  tetra- 
heda  per  minute  on  the  Cray-2.  The  computation  was  initiated  by  imposing  the  boimdary 
and  initial  conditions;  only  two  levels  of  refinement  were  used  to  refine  the  shock  at  its 
initial  position.  The  resolution  and  fidelity  of  the  simulated  shock  wave  diffraction  phe¬ 
nomena,  performed  via  the  solution  of  the  transient  compressible  Euler  equations,  were 
enhanced  by  the  application  of  the  classic  h-enrichment/coarsening  grid  adaptation  scheme, 
with  density  as  the  critical  adaptation  parameter. 

The  3-D  methodology  was  applied  in  this  project  to  the  simulation  of  sh;x:k  di&action 
processes  aroimd  several  three-dimensional  structures:  a  generic  modem  main  battlefield 
tank,  a  T-62  tank,  and  a  typical  missile  in  flight.  The  results  shown  in  this  section  demon¬ 
strate  the  successful  application  of  the  new  3-D  adaptation  procedure  to  shock  interaction 
with  curved  surfaces:  a  new  capability.  Ehccellent  shock  adaptation  and  resolution  were 
obtained  at  all  times.  All  shocks  were  captured  as  sharp  discontinuities,  without  produc¬ 
ing  pre-  or  post-shock  oscillations.  Several  interesting  three-dimensional  shock  diffraction 
processes  were  identified  and  discussed.  Among  them  are  the  time  evolution  of  three- 
dimensional  Mach-stems,  3-D  shock  wave  focusing  and  diffraction  within  comers  and 
within/between  the  wheels  and  wheels/treads/cover  assembly,  etc.  The  results  demon¬ 
strate  that  correct  load  determination  requires  the  precise  geometric  description;  simpli¬ 
fied  2-D  simulations  or  even  3-D  computations  that  use  a  simplified  geometric  description 
(such  as  representing  the  turret  eis  a  square  block)  would  yield  erroneous  results.  Among 
the  practical  results  obtained  and  discussed  were  the  increased  roll-over  probability  for  a 
side-impact  over  a  back-impact,  and  the  very  high  instantaneous  forces  exerted  on  some 
surfaces  due  to  the  high-amplitude  refiected  shocks.  Finally,  the  results  demonstrate  the 
robust  performance  of  the  method  and  show,  at  least  for  the  present  application  and  prob¬ 
ably  for  other  strongly  unsteady  fiows  as  well,  considerable  savings  in  both  CPU-time  and 
storage  over  fixed-mesh  stmctured  grid  schemes. 
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Figure  1-1.  Initial  Mesh  Refinement  Levels,  Computational  Grid  and  Pressure  Contours  for  the  Complete 
Computational  Domain  (Figures  la  to  Ic),  and  for  an  Expanded  Zone  Near  the  Box  (Figures  Id  to  If). 


Figure  1-2.  Expanded  Views  of  Pressure  and 
Den.sity  Contours  .Ground  the  Box.  t=0.4. 
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Figure  1-3.  Expanded  View  of  Pres.sure  Contours 
Around  the  Box,  t=0.6. 
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Figure  1-5.  Expanded  Views  of  Computational 
Grid,  Mesh  Refinement  Levels,  and  Pressure 
Contours  Around  the  Box,  t=0.8. 


Figure  1-6.  Expanded  Views  of  Pressure  (a),  Mach 
Number  (b)  and  Vorticity  Contours  (c)  Around  the 
Box  at  t=1.0,  and  Mach  Number  (d)  and 
Vorticity  Contours  (e)  at  t=1.2. 
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Figure  1-7.  Expanded  Views  of  Computational 
Grid  (a),  Mesh  Refinement  Levels  (b),  Pressure 
(c),  Mach  Number  (d),  and  Vorticity  Contours  (e) 
\round  the  Box,  t=1.4. 


Figure  1-8.  Expanded  Views  of  Vorticity  Con¬ 
tours  Between  the  Bottom  of  the  Box  and  Top 
of  the  Elevation;  a)  t=1.6;  b)  t=l  '  ^  t=2.0; 

d)  t=2.2;  e)  t=2.4;  f)  t=2.6;  g)  i  )  t=3.6; 

i)  t=4.2;  j)  t=5.0;  and  k)  t=6.0. 
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Figure  1-11.  Comparison  of  Experimental  and 
Numerical  Pressure  Time  Histories  and  Impnlnnf 
for  Several  Stations  Around  the  Box. 
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Figure  2-2a-d.  Expanded  Views  of  the  Mesh  Refine¬ 
ment  Levels,  Density  and  Mach  Number  Contours  at 
t=:0.125  ms,  and  Density  Contours  at  t=0.18  ms. 
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Figure  2-4a-d.  Expanded  Views  of  the  Computa¬ 
tional  Mesh,  Mesh  Refinement  Levels,  Density 
and  Mach  Number  Contours  at  t=0.35  ms. 


Figure  2-3a-e.  Expanded  Views  of  the  Density 
and  Mach  Number  Contours  at  1=0.21  ms, 
Density  Contours  at  t=0.245  ms,  and  Density 
and  Mach  Number  Contours  at  t=0.27  ms. 


Figure  2-5a.  Expanded  Views  of  the  Density- 
Contours  at  t=0.45  ms. 
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Figure  2-5b-c.  Density  Contours  at  t=0.50  ms, 
and  Close  Up  Density  Contours  at  t=0.50  ms. 


Figure  2-6a'd.  Expanded  Views  of  the  Computa¬ 
tional  Mesh,  Mesh  Refinement  Levels,  Density 
and  Mach  Number  Contours  at  t=0.60  ms. 


Figure  2-8a-e.  Expanded  Views  of  the  Mesh 
Refinement  Levels,  Computational  Mesh, 
Density,  Mach  Number  and  Vorticity  Contours 
at  t=1.90  ms. 


Figure  2-7.  Expanded  View  of  the  Density 
Contours  at  t=0.65  ms. 


Figure  2-9a-h.  Expanded  Views  of  the  FI 
Contours  at  t=0.1375  ms  and  t=0.15  ms 


Figure  2-lla-g.  Expanded  Views  of  the  Flow  Near  the  Upstream  Corner  of  the  First  Chevron;  Density, 
Mach  Number  and  Vorticity  Contours  at  t=0.45  njs,  t=:0.48  ms,  t=0.49  ms,  t=0.50  ms,  t=0.51  ms,  t=0.55 
ms,  and  t=0.60  ms. 
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Figure  3-1.  Refinement  Cases 


Figure  3-2.  Possible  Choices  for  'Inner  Diagonals’ 


Figure  3-3.  De-Refinement  Cases 


Figure  3-6a.  Expanded  View  of  Superimposed  Adapted 
Computational  Mesh  and  Pressure  Contours  on  the 
Plane  of  Symmetry,  t=4.2  ms- 
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Figure  3-7b.  Superimposed  Pressure  Contours  on  the 
Boundaries  with  Pressure  Contours  on  Two  X-Y  Planes 
(at  Z=0.5  m  and  Z=1.5  m),  t=:5.5  ms. 
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Figure  3-8a-b.  Adapted  Computational  Mesh  (a)  and 
Pressure  Contours  (b)  on  the  Tank  Surface,  t=9.67  ms. 


Figure  3-8e.  Superimposed  Pressure  Contours  on  the 
Tank  Surface  with  Pressure  Contours  on  Two  X-Y 
Planes  at  Z=1.5  m  (left)  and  Z=2.75  m  (right), 
t=9.67  ms. 


:  1 


•'■‘v 


M*‘v 


Figure  3-8f-d  Pressure  Contours  on  .\-\  Planes  Located 
at  (c)  Z=()..')  m.  Id)  Z=1.5  m  (left)  and  Z=2.(.i  m  (right). 
t=:9fi7  ms 


Figure  3-9a-b.  Adapted  Computational  Mesh  (a) 
and  Pressure  Contours  (b)  on  the  Tank  Surface. 
t  =  12.87  rns 
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Figure  S'Qc.  Pressure  Contours  on  s  X-Y  Plane, 
Z=0.5  m,  t=12.87  ms. 


Figure  3-lOc.  Pressure  Contours  on  a  X-Y  Plane 
Z=0.5  m,  t=15.53  ms. 
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Figure  3-11.  Superimposed  Pressure  Contours  on  the 
Tank  Surface  and  Pressure  Contours  on  a  Y-Z  Plane 
at  Z=2.0  m,  t=  15.53  ms. 


Figure  3-lOa-b.  Adapted  Computational  Mesh  (a)  and 
pressure  Contours  (b)  on  the  Tank  Surface,  t  15.53  ms. 


Figure  3-12a.  Pressure  Contours  on  the  Tank 
Surface,  t=19.69  ms. 


Figure  3-13a-c.  Solid  Body  Modeling  of  a 
Modern  Main  Battlefield  Tank;  Several  Views 
Around  the  Tank.  Figure  3-13d.  Surface 
Triangulation  of  the  Tank. 

Figure  3-14a-e.  Pres.sure  Contours  on  the 
Tank  Surfaces  at  200,  400,  600,  800,  and 
1000  Time  Steps,  Respectively. 
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Figure  3-14f-i.  Pressure  Contours  on  the  Tank- 
Surfaces  at  1200,  1400,  2000  and  2800  Time 
Steps,  Respectively. 

Figure  3-15a-d.  Pressure  Contours  on  Several 
X-Y  Planes:  a)  Z=0.5  m,  200  Steps;  b)  Z=1.5  m 
200  Steps;  c)  Z=1.5  m,  600  Steps;  d)  Z=1.0  m, 
800  Steps. 
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Figure  3-15e-g.  Pressure  Contours  on  Several 
X-Y  Planes:  e)  2=0. b  m,  1200  Steps;  f)  Z=1.0  m, 
1200  Steps;  g)  Z=1.5  m,  1200  Steps. 

Figure  3-16a-d.  Superposition  of  Pressure 
Contours  and  Adapted  Grids  on  the  Tank 
Surface,  Back  Impact,  at  200,  400,  600  and 
1000  Time  Steps,  Respectively. 
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Figure  3-17a-c.  Superposition  of  Pressure 
Contours  and  Adapted  Grids  on  the  Tank 
Surface,  Side  Impact,  at  0,  300,  and  600 
Time  Steps,  Respectively. 

Figure  3-18a-d.  Pressure  Contours  on  the 
Tank  Surfaces:  a)  Tank,  300  Time  Steps; 
b)  Expanded  View  Between  the  Wheels, 
300  Time  Steps;  c)  Tank,  600  Time  Steps; 
d)  Tank,  800  Time  Steps. 
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OVERPRESSURE,  NON-DIMENSIONAL 


0123456  789  10 


TIME,  non-dimensional 


TIME,  NON-DIMENSIONAL 

Figure  3-19^,  Comparuon  of  Presrare-time 
Histories  at  Several  Locations  on  the  Tank; 

Back  Impact  Vs.  Side  Impact;  a)  Station  64, 
Back-end,  Plane  of  Symmetry;  b)  Station  29, 
Back  I^k,  Plane  of  Symmetry,  Close  to  Turret; 
c)  Station  34,  Back  Deck,  Plane  of  Symmetry, 
Further  Upstream;  d)  Station  43,  Back  Deck; 

e)  Station  57,  Back  Deck,  Further  Upstream; 

f)  Station  25,  on  T\irret;  g)  Station  10, 

'Dirret,  FVont;  h)  Station  3,  FVont  Deck; 

i)  Station  84,  Bottom,  Plane  of  Symmetry. 
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